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COHORT-II - A MONTE CARLO GENERAL PURPOSE SHIELDING COMPUTER CODE 
by Leonard Soffer and Lester Clemons, Jr. 

Lewis Research Center 

SUMMARY 


COHORT -II is a Monte Carlo general purpose shielding computer program written 
completely in FORTRAN -IV for the IBM-7094 computer. In its present form, it repre- 
sents a major revision and modernization of the original COHORT code. It has been in- 
tegrated into a single package with more flexible geometry and biasing. 

Given an initial source distribution of neutrons or photons, the code calculates num- 
ber fluxes for these particles in energy bins for arbitrarily located point detectors and 
track length fluxes averaged over the regions of a reactor -shield geometry as specified 
by the user. Energy integrated quantities such as dose rates are also obtainable. 

The code method employed is the conventional one where individual particles are 
generated, tracked, and analyzed. The shielding geometry allowed is flexible with space 
divided into homogeneous regions bounded by surfaces describable by a general quadra- 
tic equation. Spherical, cylindrical, and rectangular parallelepipeds are allowed for the 
source regions. Biasing options include selection of source particles from preferred 
locations, energies, and directions. The exponential transform can also be used during 
particle tracking. Secondary gammas can be generated and placed on tape during a 
neutron problem to be tracked during a subsequent computer run. Variable dimension- 
ing is employed for efficient use of computer storage. 


INTRODUCTION 


COHORT-n is a Monte Carlo general purpose shielding computer code. Given an ini- 
tial source distribution of neutrons or photons, it computs number fluxes for these source 
particles at arbitrarily located points or track length fluxes averaged over the regions of 
a reactor-shield geometry. It is written in FORTRAN-IV for the IBM- 7094 computer. 

In its present form, it represents a major revision and modernization of the original 
COHORT code (refs. 1 and 2). The code CAVEAT (ref. 3) similarly has been derived 
from COHORT. This report is intended as a guide to assist potential users in data pre- 


paration and explanation of the COHORT -IE code method and is intended to stand alone 
without requiring the user to refer to references 1 and 2. Those interested in CAVEAT 
should consult reference 3. 

Many features of the original code have been retained such as the specialized input 
libraries which are used to describe the geometry, the cross sections, and the use of the 
exponential transform. In addition, the collision mechanics (except for low energy neutron 
scatter) and the basic scoring methods in the original code, including statistical estima- 
tion of fluxes at point detectors and track lengths in volume detectors, have been retained 
largely intact. 

The principal modifications and additions made in COHORT -II include (1) an exten- 
sion of the geometry capability to allow for reactor -shield geometries that may be de- 
scribed by a general quadratic surface (with suitable choice of the coefficients), (2) ad- 
dition of source particle biasing to allow preferential selection of source particle co- 
ordinates, directions, or energies, (3) introduction of an estimate of the statistical 
precision for the calculated results, (4) integration of the separate segments of the code 
to allow a complete problem, from source particle generation to output of fluxes at point 
detectors and in regions, to be run in a single pass on the computer, and (5) addition of 
variable dimensioning to allow for more efficient use of limited computer core storage. 

A brief description of the COHORT-H code follows. 

The code can track either primary neutrons or photons born in rectangular paral- 
lelepiped, spherical, or cylindrical source regions that are oriented parallel to the co- 
ordinate axes. Primary particles are generated from previously known spatial, angular, 
and energy distributions supplied by the user. Secondary photons born anywhere within 
the reactor -shield geometry can also be tracked in a subsequent run. The shield geome- 
try is more flexible than the primary source geometry. Space is divided into homogene - 
ous regions bounded by surfaces that can be described by a general quadratic equation 
with suitable choice of coefficients. This allows a description of such surfaces as cylin- 
ders, cones, ellipsoids, planes, etc. , at an arbitrary orientation to the coordinate axes. 

The code requires input cross-section data in its own specialized format and can- 
not accept multigroup cross-section data as used in conventional transport codes. Pre- 
paration of cross-section data, particularly for neutrons, is tedious because it has 
not yet been completely automated. A library of neutron cross sections in the required 
format has been compiled for about a dozen elements and is available with the code. De- 
tails of this library are given in the section CROSS SECTIONS. Preparation of photon 
cross-section data is considerably simpler and can easily be done by the user. 

Biasing features of the code include preferential source particle selection and the 
exponential transform, or path -stretching, during particle tracking. The degree to 
which the exponential transform is applied can be varied depending on the direction of 
the particle and the region it is located in. 
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Particles are tracked one at a time. Histories are terminated when a particle 
escapes from the system under consideration, when its energy or weight falls below a 
previously assigned cutoff value, or when it has undergone more than a previously 
specified number of collisions. When secondary gammas are to be generated, neutron 
collision data are written on a tape and stored for future use. If leakage particles are 
to be used as a source in a subsequent problem, a leakage tape can be generated for this 
purpose. 

Code printed output includes number fluxes that are averaged over regions and 
fluxes that are at arbitrarily located point detectors. The volume averaged region fluxes 
can be sorted according to arbitrary energy bins. The point -detector fluxes list uncol- 
lided (for isotropic sources only) and scattered flux contributions separately and can be 
sorted according to their polar and azimuthal angles of arrival at the detector as well 
as their energies. An estimate of the statistical precision is included for both volume 
and point detectors. Energy integrated quantities such as dose rates, heating rates, re- 
action rates, and fluxes greater than a given energy can also be computed for both volume 
and point detectors by entering suitable response functions. Where only one type of de- 
tector is of interest, the computation associated with the other type can be omitted to 
save time. All scoring, both for volume and point detectors is performed by statistical 
estimation. 

The program does not require the user to handle tapes except when a leakage or a 
secondary gamma tape has been requested. A complete analysis of a typical reactor - 
shield problem would require three separate Monte Carlo problems to be run on the 
computer. These runs would be used to track and analyze, respectively, primary neu- 
trons, primary photons, and secondary photons. A physical source tape for the second- 
ary gamma problem would be generated during the primary neutron run. 

The sections that follow describe the code method and input data requirements in 
detail. Sample problem input and output are discussed in appendix A. Flow charts for 
each subroutine in the code are also included and shown in figures 5 through 23. 


CODE METHOD 
GEOMETRY 

The geometry of a problem is defined by identifying regions and their bounding sur- 
faces. A region is a localized part of three-dimensional space bounded by one or more 
surfaces (up to a maximum of 12 boundaries per region). The nuclear composition of 
the material within a region must be the same for every point within it. 

The geometry calculation performed by the code has two basic functions. Given the 
spatial coordinates and direction cosines of a particle, these are (1) to determine the 


3 



region the particle is located in, and (2) to compute the distance from the particle to the 
nearest boundary in its direction of travel. 


Geometry Method 

The geometry method used is the same as that of the QAD point -kernal shielding 
computer code (ref. 4) and will be described only briefly. The shapes of the regions 
that can be described are limited to those whose bounding surfaces can be expressed as 
special cases of a general quadratic equation. In the code the basic equation and four 
special forms are included to simplify input data. These equations are as follows: 


2 2 2 

Ax + By + Cz + Exy + Fxz + Gyz + Hx + Jy + Kz - D = r (1) 

A (* - x Q ) 2 + B(y - y Q ) 2 + C(z - z Q ) 2 - D = r (2) 

x - D = r (3) 

y - D = r (4) 

z - D = r (5) 


Equation (1) is the general quadratic equation and can be used to describe such con- 
vex surfaces as spheres, cylinders, cones, ellipsoids, hyperboloids, and planes at ar- 
bitrary locations and orientations to the coordinate axes by suitable choices of the coef- 
ficients. Equation (2) can be used to describe spheres, cylinders, ellipsoids, and hyper- 
boloids parallel to the coordinate axes. Finally, equations (3) to (5) represent planes 
perpendicular to the x, y, and z axes, respectively. The user should try to select the 
simplest equation needed to describe the required surface since this not only reduces 
required input data, but computer time as well. 

The variable r, also known as the residual, defines how an arbitrary point in space 
with coordinates (x, y, z) is related to the surface. The expression r(x, y, z) < 0 defines 
all points on one side of the surface, and r(x, y, z) > 0 defines those points on the other 
side. The residual is zero for all points on the surface. For each surface bounding a 
region, a direction index j that has a value of either +1 or -1 must be defined. The 
value of the direction index is chosen so that the product of j and the residual r(x, y, z) 
is negative if the point (x, y, z) lies inside the region. This is the same convention used 
in the QAD code. Each j must be chosen properly according to the type of bounding 
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surface. If the bounding surface is closed (i. e. , sphere, ellipsoid, etc. ), a region in- 
side the surface would have j = +1, and a region outside would have j = -1. An open 
bounding surface (i. e. , cylinder, paraboloid, cone, etc.), would have j = +1 for a re- 
gion between the surface and the axis of the figure and j = -1 for all other space. A 
plane surface would have j = +1 for a region to the left side of the plane and j = -1 for 
a region to the right side. 

A single bounding surface may be used in the definition of many regions. The re- 
gion definition requires the number of bounding surfaces (NB) of that region. If NB is 
entered with a positive sign, this indicates that the region is an inside region, while a 
negative sign affixed to NB indicates that the region is an outside region. Once a parti- 
cle enters an outside region its history is terminated. Thus, regions that are of no in- 
terest in particle tracking can be designated as outside regions and their outermost 
boundaries need not be defined. An outside region should not be embedded within one or 
more inside regions, however, because particles will not be tracked across the outside 
region. 

The region definition also requires a list of each bounding surface (IBN) that the 
region comprises. The sign affixed to the bounding surface IBN serves as the j index 
discussed previously. In addition, for each IBN listed, the user must also read in the 
most probable region a particle will enter when it leaves the region being defined. 

To illustrate this, consider the simple example of describing a right circular cylin- 
der bounded by two parallel planes. A sketch is shown in figure 1. The regions are in- 
dicated by numbers enclosed in a circle, and the boundaries are indicated by numbers 
enclosed in a square, one line of which is the boundary involved. A sketch of the prob- 
lem geometry, with suitable identification of regions and boundaries, is strongly recom- 
mended as the starting point for making out the geometry input. 

Boundaries 1 and 2 are planes parallel to the z-axis (type 5), located at z =0 and 
10, respectively; boundary 3 is a cylinder parallel to the z-axis (type 2) with a radius of 
5. Region 1 consists of all points that have z values (0 s z < 10) and are simultane- 
ously within the cylinder. The points of region 2 have the same z values but are out- 
side the cylinder. Regions 3 and 4 are the sets of points with z < 0, and z > 10, 
respectively. 

We now wish to determine the direction indices of region 1 for each of its three 
boundaries. The direction index of boundary 1 with respect to region 1 is -1. This is 
because boundary 1 has the form z - 0 = r. All points in region 1 have positive z, 
thus the residual r for boundary 1 is positive. Since the product of the residual and 
direction index must be negative, this is achieved by setting j = -1. By the same 
reasoning, the direction index of boundary 2 with respect to region 1 is +1, and the 
direction index of boundary 3 is +1. 
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In defining region 2, it can be seen that it has the same three boundaries as re- 
gion 1. The direction indices of boundaries 1 and 2 with respect to region 2 have the 
same values as those with respect to region 1. But since region 2 defines those points 
outside the cylinder, the direction index of boundary 3 with respect to region 2 is now 
-1. Regions 3 and 4 are easily defined. Each has but one boundary. The direction 
index of boundary 1 with respect to region 3 is +1, and the index of boundary 2 with 
respect to region 4 is -1. For detailed input instructions, see the section Input Data. 


Region Location 

To determine the region a particle is located in, the code commences searching 
through the regions, one at a time. For each region, the product of the residual r and 
the direction index j is computed for every one of its boundaries. If this quantity is 
negative for every boundary in the region, the particle is concluded to be in the region. 

If the product is positive for any boundary of a region, the search in that region is 
abandoned immediately and proceeds to the numerically next higher region. If the code 
cycles completely through the regions and fails to place the particle in one of them, a 
geometry error is printed out and the problem is terminated. When a particle crosses a 
boundary from one region to another, the region search begins with the most probable 
region entered as defined in the region description. 


Distance to Nearest Boundary 

Once a particle has been located in a region, the distance to the nearest boundary 
of the region along its direction of travel is then computed by checking the distance to 
all the boundaries of the region. If the particle coordinates are x’, y', and z', if its 
direction cosines are a, (S, and y, and if the distance to the nearest boundary along its 
direction of travel is S, the coordinates of intersection (x, y, z) of the particle’s path 
with the boundary are 


x = x* + aS 

y = y’ + /3S (6) 

Z = Z’ + yS 
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Since the point of intersection x, y, z is simultaneously on the particle flight path 
and on the boundary, equations (6) and the appropriate boundary equation can be solved 
simultaneously to yield an equation which in general is quadratic in S. This equation 
can be written in the form 


P 1 S 2 + 2P 2 S + P 3 a 0 (7) 

where P^, P 2 , and Pg are functions given in terms of the boundary equation coeffi- 
. cients and particle direction cosines. This equation is then solved for S. The physical 
intersection of the particle with the boundary is given by the smallest, nonnegative, real 
value of S. Then S is calculated for all the other boundaries of the region and com- 
pared. The smallest, nonnegative, real value is retained as the distance to the nearest 
boundary in the direction of travel. A small quantity EPSL is added to step the particle 
across the boundary. If a particle is located in a given region, but the routine is unable 
to calculate a satisfactory distance to one of its boundaries, an error message is printed 
out and the problem is terminated. 


CROSS SECTIONS 

This section is intended to describe the type of cross-section data required and how 
the data are handled in the code. As mentioned in the INTRODUCTION, a library of 
neutron cross sections in the required format has also been prepared and is available 
with the code. The elements in this library are hydrogen, lithium 6 and 7, carbon, 
oxygen, sodium, aluminum, chlorine, calcium, iron, tungsten, and lead. The sources 
of these data are noted on the comments cards included with the code. Additional details 
about this library are given in the next section. 


Input Requirements 

Cross-section input data for tracking neutrons are entered into the code in three 
groups of data designated as library numbers 6, 9, and 10. For gammas, only library 6 
data are required. In both cases, these are entered as microscopic cross sections, with 
material or atom densities entered separately. Macroscopic cross sections are com- 
puted internally. 

Library 6 data are the various reaction cross sections at discrete energies with 
energy points spaced arbitrarily by the person generating this table. Actual values of 
the cross sections for a particle at a specific energy are obtained by linear interpolation 
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within the table. The interpolation routine performs efficiently by progressively halving 
the energy range of the search. The elements in the library 6 neutron data package, 
available with the code, cover the energy range from 18 MeV to 0. 025 eV in 90 discrete 
energy points, with a resolution of 10 energy points per decade. This resolution is be- 
lieved to be adequate for those shield problems where an integral answer such as total 
dose rate or heating is of primary concern; but is inadequate for those problems where 
highly detailed spectral information is desired. 

For neutrons this table lists the energy, the total cross section, the elastic plus 
inelastic scattering cross section, and the elastic scattering cross section. Those 
reactions resulting in the disappearance of a neutron such as (n, y), (n, a), (n, p), etc. , 
are represented by internal computation of the total cross section minus the sum of the 
elastic and inelastic scattering cross section. Neutron interactions considered are 
neutron absorbtion (or disappearance), elastic scattering, and inelastic scattering. The 
(n, 2n) reaction is handled in the data package as an inelastic scattering. 

’ *“*For gammas, library 6 lists the energy, total cross section, Compton Shattering 
plus pair -production cross section, and the Compton scattering cross section. Gamma 
interactions considered are photoelectric absorption, pair production, and Compton 
scattering. Photoelectric absorption of gammas is computed as the difference between 
the total and the sum of the Compton plus pair -production cross sections. No library 
of gamma cross sections has been included with the code, because of the relative ease 
of preparation of these data from available compilations. However, a gamma cross- 
section library is available on request. 

Library 9 data are required for those elements undergoing neutron inelastic scatter- 
ing. The energy levels E. of the nuclide, a list of incident neutron energies E-, and a 

J ^ 

two-dimensional array P^ are required for each element. The array P^ represents 
the cumulative probability of a neutron of incident energy exciting the nuclide to an 
energy level Ej. The probability for neutrons at an arbitrary incident energy E is 
obtained from the probabilities for the largest value of E^ less than E. When an in- 
elastic scatter occurs, the energy level Ej to which the nuclide is excited is computed 
by choosing a random number and finding the smallest value of j that satisfies 

E p ii a R i 

j 

Library 10 data are required for all elements except hydrogen. This library de- 
scribes the anisotropic neutron elastic scattering in the center -of -mass system. A list 
of incident neutron energies is entered. A two-dimensional array listing cosines of the 
scattering angle (center of mass) for each incident neutron energy is also required. The 
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values of these cosines are obtained by integrating the differential elastic scattering 
cross section and finding values of the cosines of the scattering angle that divide the in- 
tegral into a number of equal parts. The cosine values that define the histogram are the 
bulk of the library 10 data. The number of segments of the histogram can be varied to 
adequately represent the differential scattering cross section. The library data included 
with the code divide the integral into 20 equal parts. When an elastic scattering has 
taken place and the polar angle of scattering is required, a random number is chosen. 
Because each segment of the histogram represents a probability Pj of the neutron hav- 
ing scattered through that angle, the segment where scattering occurs for a given history 
is computed as 


i> i*h 

i=l 

Once the segment is defined, a second random number is chosen, and the exact angle of 
scattering is located at random within the segment. 

The remaining required input data to describe the nuclear properties of the materials 
are the atom densities and the atomic weights (not needed for gammas). 


Internal Data Preparation 

After the cross sections have been read in, the macroscopic material cross sections 
are computed for neutrons or gammas using the atom densities and the library 6 tabula- 
tions. Certain probabilities are then computed for use in later calculations. For neu- 
trons, the probability of a scatter collision being an elastic scattering event P g ^ is de- 
fined to be 


p _ CT el 

F el 

°el + °inel 

and is computed for each element. For gammas, similarly, the probability that a given 
scatter is a Compton scatter 


g compt 

CT compt + °pair 
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is computed. It should be noted in this regard that pair production can be regarded as a 
type of scatter since photons resulting from the pair annihilation emerge from this event. 


Cross-Section Requirements for Secondary Gammas 

Input cross sections for capture gamma source generation consist of a set of tables 
for each element where capture gammas are possible. These tables include (1) a list of 
incident neutron energies, (2) a list of gamma energies, and (3) a two-dimensional array 
Nj^, giving the number of photons of gamma energy E^. produced by an incident neutron 
of energy Ej. 

Since inelastic gammas are determined from the library 9 data, no additional cross 
sections are required to generate the inelastic gamma ray sources. 


PRIMARY SOURCE PARTICLE GENERATOR 
General 

The primary source particle generator is a procedure for generating nine parame- 
ters used to define each primary neutron or photon. These parameters are the three 
spatial coordinates of its birthplace, three direction cosines of its initial direction, its 
initial energy, weight (different from unity if source biasing is performed), and region 
of birth. The parameters are written on a tape consisting of 100 particles per record 
and are used as input to the particle tracking routines. 

The basic assumption of the method is that the total distribution defining the source 
can be separated into distributions which define the spatial, direction and energy para- 
meters independently of one another. In addition, the spatial distribution is assumed to 
be separable and can be defined as the product of three independent, single- coordinate 
distributions. 

In keeping with these assumptions, the code link SOURCE contains four subroutines 
labelled SOURCE, DIRCOS, COORD, and ENERGY. SORCE is the main routine and uses 
the others to independently select, respectively, direction cosines, spatial coordinates, 
and the energy of each particle. These are then combined and a suitable weight is com- 
puted. 

Source geometry options available to the user include specification of a point source 
or of a rectangular parallelepiped, cylindrical, or spherical volume sources. Angular 
options include monodirectional sources or isotropic and anisotropic sources. Energy 
options include monoenergetic or arbitrary energy spectrum sources. 
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A useful tool available to the user is the ability to perform biasing in the selection of 
any or all of the source particle parameters. For example, source energy biasing may 
be accomplished by reading in two energy distributions: the actual one and a biased 
distribution. The particle energy is selected from the biased distribution, but its weight 
is computed as the ratio of the actual to the biased distribution, thereby removing the 
bias. Similarly, directional source biasing can be performed, for example, for an iso- 
tropic angular distribution, by starting more particles, each of a lower weight, in a pre- 
ferred direction than would normally come from an unbiased isotropic source. Biasing 
the initial polar (measured from the z-axis) or azimuthal angle can be performed inde- 
pendently, if desired. Finally, the selection of the particles initial spatial coordinates 
can also be biased. Biasing for one coordinate can be performed independently of the 
others. 


Source Spatial Coordinate Selection 

The initial spatial distribution of the primary particles, indicated by the control 
parameter IC, may be for a point source or for rectangular parallelepiped, cylindrical, 
or spherical volume sources (with or without biasing). 

The option IC = 0 indicates that a point source is desired. Three spatial coordi- 
nates x, y, and z are entered and all source particles are born at this point. A spatial 
weight W g p of 1 . 0 is assigned to each particle. 

The options IC = 1, 2, or 3 refer, respectively, to the rectangular parallelepiped, 
cylindrical, and spherical volume source geometries. The spatial distributions and 
coordinate values are entered, in pointwise fashion, for each of the three coordinates 
defining the source. For the users convenience, the input spatial distributions need not 
be the actual probability density functions, but may be relative distribution functions. 

The code reads in all the relative spatial distributions, integrates them over their re- 
spective coordinates by trapezoidal integration, and normalizes them internally to obtain 
actual probability density functions. 

Having obtained the probability density function f (77) for each coordinate, the proba- 
bility Pp that a particle’s coordinate n is located in the v" interval between 77 . and 
7?i + l is given by 



f(?7) d-q = P A 
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A random number Rj is chosen, and the coordinate interval is selected by talcing the 
smallest value of i that satisfies 


Z p i £K i 

i 

The actual value of the coordinate r\ is determined by selecting a second random num- 
ber and evaluating 77 at random within the interval. A spatial weight of 1 . 0 is also as- 
signed to each particle. 

When spatial biasing is to be performed, a negative sign is affixed to the IC para- 
meter. In addition to the actual distributions for each coordinate, the user must also 
read in biased distributions for all coordinates, even though spatial biasing is to be 
performed for only one. If the biased and the actual spatial distributions are identical 
for a given coordinate, then, effectively, no biasing in the selection of that coordinate is 
performed. As with the actual distributions, the biased distributions entered may be 
relative ones and are normalized internally to obtain biased probability density functions 
f *('/]) . The coordinate interval and value are selected in a similar manner as that dis- 
cussed previously, except that samples are taken from the biased probability density 
function P*. A coordinate weight W c is computed as 


/Vi 

J Hn) &n 



and a spatial weight W is computed as the product of the three coordinate weights. 

sp 


Source Direction Cosine Selection 

The initial angular distribution of the primary particles, indicated by the control 
parameter, ID, may be monodirectional, isotropic (with or without biasing), or aniso- 
tropic. 

The option ID = 1, indicates that the source particles are monodirectional. A set 
of direction cosines a, /3 , and y are input to the program, and all source particles 
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are emitted with these specified direction cosines. An angular weight of 1. 0 is assigned 
to each particle. 

The option ID = 2 is used when source particles are emitted isotropically. The 
direction cosines are determined by a rejection technique that was developed by Kahn 
(ref. 5) and is shown in figure 2. Each particle is given an angular weight of 1. 0. No 
additional input is required for this option. 

If particles are to be emitted isotropically but biased in either the polar or azimu- 
thal angular directions, the option ID = -2 is used. The cosines defining each polar 
angular interval and the probability of particle emission within each interval are read 
in. Similarly, the angles defining each azimuthal angular bin and the probabilities of 
particle emission in each is also required. The cosine of the polar angle is selected by 
obtaining a random number Rj and choosing the smallest value of i for which 

^P. (cos e) ^ Rj 
i 

where P^(cos 0) is the probability (biased) of particle emission in the i n polar angular 
interval. 

With the polar angular interval known, the cosine of the polar angle y is chosen 
at random within the interval by requiring that 

y = ri + R 2 ( n+i ' 

■fVi 

where y ^ and y. + ^ are the values of the cosines defining the i n interval and Rg is a 
random number. 

A polar angular weight is then calculated as the unbiased divided by the biased 
probability of emission per interval as 


W - ^ ~ y l+l 
9 2P i (cos e) 

Calculation of the azimuthal angle is done similarly by choosing the smallest value 
of i such that 



sR 3 
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The azimuthal angle is then computed as 


cp= <p i + R 4 (<p i+1 - <p$ 
An azimuthal angular weight is computed as 



gtfi - n 

2jrPj(«p) 


An angular weight is computed as the product of the polar and the azimuthal weights. 
The direction cosines a and /3 are given by 


and 




cos cp 




sin cp 


The option ID = 3 is used to specify an anisotropic angular distribution. Input is the 
same as that for the ID = -2 option. The selection of polar and azimuthal angles are 
performed identically, but an angular weight of 1. 0 is assigned to each particle. 


Source Energy Selection 

The initial energy distribution of the primary particles, indicated by the control 
parameter, IE, may be monoenergetic or may be selected from an arbitrary distribution. 

The option IE = 1 indicates that the particle energies are to be monoenergetic and 
the energy value (E) is entered as input. An energy weight of 1. 0 is assigned to each 
particle. 

Option IE = 2 is used when energies are selected from a continuous energy spectrum. 
The spectrum is divided into intervals. The probability that a particle is emitted in the 
i m interval P^(E) is read in for every interval. The interval of the source particle is 
selected by choosing a random number and computing the smallest value of i for 
which 
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E p i (E) a R i 

i 

The actual particle energy is obtained from 


E = Ej + R 2 (E i+l " E i } 

where Ej is the lower bound of the i 1 energy interval. An energy weight of 1. 0 is 
assigned. 

Option IE = -2 allows the user to perform energy source biasing. As before, the 
probability of particle emission in the i energy interval P^(E) is required for every 
interval. In addition, the biased probability emission per interval P*(E) is also re- 
quired. The code determines the energy of the source particle the same as in option 
IE = 2, except that sampling is performed from the biased probability distribution. An 
energy weight is calculated as 


W 


E “ 


P t (E) 

P*(E) 


Stratified Selection of Source Particle Parameters 

Since the distribution of the primary source particles in space, direction, and energy 
is known, a stratification procedure can be used to assign the source particle parame- 
ters. If the probability of a particle being located in interval i is P., then the total 
number of particles to be located in the i m interval is NPp where N is the total num- 
ber of source particles to be generated. This procedure, known as stratification, is 
more desirable than choosing the interval by comparing a random number with the cu- 
mulative probability distribution. The stratification procedure will always select ex- 
actly NPj particles within the interval; whereas, the probability selection procedure 
only does so on the average. 

Stratification in the selection of all of the primary source particle parameters is 
impractical, however. Imagine a typical problem where the spatial coordinates are 
divided into 10 intervals for each coordinate, the angular coordinates are divided into five 
intervals each, and the energy spectrum is divided into 20 intervals. If the selection pro- 
bability per interval is equally likely, the probability of a particle being located in a cell 
defined by the coordinates x, y, z, 9, cp, and E becomes 


15 



To locate just one particle per cell would require generation of about 10 source parti- 
cles. If the probabilities for each cell are not equal, as is more likely, assuring at 

least one particle per cell requires even larger N. Monte Carlo calculations are time 

4 5 

consuming, and typical problems involve generation of about 10 to 10 primary parti- 
cles. 

By proper designation of an input parameter (see the section Input Requirements for 
details), the code allows any one of the primary particle source coordinates such as 
energy, polar angle, etc. , to be selected for stratification. 

The number of particles located within an interval of a stratified coordinate is NP., 
where is the probability of a particle being born within the interval. The intervals 
for the remaining coordinates are selected by comparing a random number with the 
cumulative probability distribution of that coordinate. No more than one coordinate can 
be stratified, or, if the user wishes, none of the coordinates need be selected by strati- 
fication. 


PARTICLE TRACKING AND SCORING 

The particle tracking and scoring routines compose the third link of the COHORT -II 
program. For convenience, the discussion is divided into particle tracking and particle 
scoring. In the code, a given routine may perform part of, both functions, for greater 
overall efficiency. Four routines are in this link. They are HTA, HISTRC, SCATR, and 
GEOMT. HTA acts as an overall executive routine in the particle tracking, but it also 
computes the uncollided fluxes at point detectors and performs preliminary statistical 
tallying. HISTRC is concerned primarily with particle tracking between collisions, and 
it performs the exponential transform and calculates fluxes within regions. SCATR 
handles the collision processes, energy, and directions after scatter and also computes 
scattered fluxes at point detectors. Finally, GEOMT is concerned only with geometry 
calculations and has been discussed previously. (See GEOMETRY section. ) 

The following sections discuss these functions in more detail. 


Particle Tracking 

Distance to collision . - The number of mean free paths D2 T (E) that a particle will 
travel on its flight to its next collision is selected from the exponential distribution 



P [x < DZ t (E) < (x + dx) = e' x dx 


by taking 

-In R n 

D = - 

s t (e) 

where is a random number, and S,p(E) is the total macroscopic cross section. 

After D is computed, a check is made to determine whether the flight path is con- 
fined to one region, in which case calculation of the collision location is straightforward. 

If the particle leaves the region, the exit position is determined and the path length S 
in the region is computed. 

The distance the particles can travel in the new region entered is computed by the 
equation 

£ t (E) 

D’ = (D - S) -f— 

z T (E) 

where S^(E) and S,p(E) are the total macroscopic cross sections for the region entered 
and the region just crossed, respectively. If D f is less than the distance across the new 
region, the collision occurs within it. If D* is greater, the process is repeated as 
often as necessary. 

Coll ision process. - This discussion covers how the types of collisions undergone 
by neutrons and gamma rays are treated in the code and the determination of the energy 
and direction after scattering. 

Only scattering collisions are treated specifically in the code. These include elastic 
and inelastic scattering for neutrons and Compton scatter and pair production for gamma 
rays. To account for particle absorption, the weight of the particle after a collision is 
computed as the product of the weight before collision and the probability that the colli- 
sion was a scattering event. Thus, 

W f = WP d (E) 

b 

where W and W* are the particle weights before and after collision, respectively, and 
PtE) is the scattering probability for particles of incident energy E. The sequence of 

o 

events in the code in handling the collision processes for both neutrons and photons is 
as follows: (1) the scattering element or nuclide is determined, (2) the particle is tested 
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to determine whether it is a neutron or a photon, (3) the type of scatter event is deter- 
mined (elastic against inelastic for neutrons; Compton against pair -production for pho- 
tons), (4) the polar scattering angle and the energy after collision are found using the 
appropriate scattering laws, depending on the type of particle and type of scatter event, 
(5) the azimuthal scattering angle is computed, and, finally, (6) new direction cosines 
after scattering are computed from a knowledge of the values before scattering and the 
polar and azimuthal scattering angles. These six steps are discussed in further detail 
in the remainder of this section. 

When a collision occurs, first the scattering nuclide is selected by determining the 
smallest value of i for which 


1 

S S (E) 


£ Vsi< E > 2 R 


where £ C (E) is the macroscopic scattering cross section of the material, a .(E) is the 
o th S1 

microscopic scattering cross section for the i element, and N. is the atom density 

for the l element. A test is then made to determine whether the particle is a neutron 

or photon. 

For gamma rays, the type of scattering event is determined by comparing a random 
number to P (E), the probability that a given scatter is a Compton scatter. If 
P c (E) 5: R^, the event is a Compton scatter. Otherwise, pair production results. If the 
event is a Compton scatter, the energy and direction after scattering are selected from 
the Klein -Nishina formula by use of a rejection technique developed by Kahn (ref. 5). 
This rejection technique is illustrated in figure 3. If a pair production event occurs, it is 
assumed that the electron -positron pair are annihilated immediately, resulting in two 
isotropically emitted photons, each of energy 0. 511 MeV. Only one photon is followed 
in the code, but its weight is doubled. The direction cosines of the isotropically emitted 
photon are obtained via a standard rejection technique. 

For neutrons, the type of scattering event (elastic or inelastic) is chosen by picking 
a random number Rj and comparing it to P gl (E), the probability that a given scatter 
is elastic. If P g j(E) ^ Rp the event is taken as an elastic scatter. If not, the scatter 
is inelastic. 

The polar scattering angle in the center -of -mass (CM) system is then determined. 
Neutron elastic scattering is assumed to be isotropic in the center -of -mass system if 
the neutron incident energy is less than an input parameter ESI(J), which is entered for 
each of the J elements. If the incident energy is greater than ESI(J), elastic scattering 
is taken to be anisotropic. For isotropic elastic scattering, the polar scattering angle 
in the center -of -mass system X is computed from the equation 
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cos X = 2R - 1 


where R is a random number. For anisotropic elastic scattering, the cosine of the 
polar scattering angle in the center -of -mass system is determined from input cumula- 
tive probability distribution tables (library 10 data) as a function of incident neutron 
energy and element. 

Inelastic neutron scattering is assumed to be isotropic in the center -of -mass system 
for all elements. The polar scattering angle is computed as 

cos X = 2R - 1 


where R is a random number. 

The energy after scattering is then calculated. A parameter B is assigned a value, 
depending on the type of scattering event. For neutron elastic scattering, B = 1. For 
neutron inelastic scattering, B is defined as 


B = 



'A+l\ E x 
, A / E 


where E is the excitation energy which is selected from a probability distribution 
x 

table (library 9 data) for each element, E is the neutron incident energy, and A is the 
atomic weight for the element. The energy after scattering E* is then computed to be 

E’ = — (l + A 2 B + 2AB 1 / 2 cos x) 

(A + l) 2 


where the same equation is used for both elastic and inelastic scatter events, substituting 
the proper value for B. Use of this equation with different definitions of B, depending 
on elastic or inelastic scatter, is taken from the original versions of the code (refs. 1 
and 2) and has been retained for economy of programming and convenience. 

The polar scattering angle in the laboratory system i // is then calculated from the 
center -of -mass system angle by means of the equation 


cos i// = 


(1 + AB^ 2 cos X) 


(l + A 2 B + 2AB 1//2 cos x) 


1/2 
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where again the same equation is used for both elastic and inelastic scattering by sub- 
stituting the appropriate value for B. 

A simple case arises if the element is hydrogen. Neutron scattering in hydrogen is 
assumed to be isotropic in the center -of -mass system, and in addition there is no in- 
elastic scattering. The preceding equations take on a particularly simple form since 
A = B = 1. The equation for the polar scattering angle becomes; 

cos X = 2R - 1 


for the energy after collision, 


E’ = - (1 + cos X) 
2 


and, for the laboratory scattering angle, 


cos \j/ = 


/ 1 + cos X 

V 2 ~ 


1/2 


Eliminating (1 + cos X) from the second and third equations and substituting the first 
equation into the third, we obtain 


2 

E* = E cos 

and 

cos i// = Vi 

The cosine of the laboratory scattering angle is computed somewhat more quickly by 
restoring to a scheme developed by Kahn (ref. 5). Two random numbers are chosen, 
and the largest of these is taken to be equal to the cosine of the scattering angle. Thus, 

cos = Vi = max (R^> R 2 ) 

Another special case arises if the neutron energy is below an input value ETM, 
where ETM is assumed to be a thermal energy level. Neutrons with energy less than 
ETM are scattered isotropically in the laboratory system without energy change. 
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For those particles not isotropically scattered in the laboratory system (Compton 

scattered photons and neutrons with energy greater than ETM), the azimuthal scattering 

angle cp is next selected from a uniform distribution in the interval from 0 to 2n. A 

proper choice of cp would be 27 tR, where R is a random number. However, since 

sin cp and cos cp are desired, a commonly used rejection technique is employed for ef- 

2 2 

ficiency. Two random numbers are chosen, and the sum of the squares S = R^ + Rg is 
tested. If S < 1, then 


cos cp = 



and 


sin cp = 


2RjRg 


If S > 1, two other random numbers are chosen and the cycle is repeated. 

The direction cosines of the particle after collision is computed from a knowledge of 

the polar and azimuthal scattering angles ip and cp in the laboratory system and from 

the direction cosines a, p, and y of the particle before collision. The new direction 

o 

cosines are computed as follows: when (1 - y ) > 0. 0001, 

a , = ay sin xpcos cp - p sin ^ sin cp + a cos ^ 

(1 - r 2 ) 172 

s t = a.sin> sin cp + y p sin x}/ cos cp cog 

(l -r 2 ) V2 


and 


1/2 


/ 2 \ ' 

y 1 = y cos ip - II - y ) sin ip cos 




when (1 - / ) < 0. 0001, 


o?’ = sin \p cos cp 


P* = sin ip sin cp 
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and 


y* = cos i// 

Exponential transform . - The exponential transform finds its greatest application in 
deep penetration problems, where it effectively serves to stretch the path length of the 
particle between collisions. To achieve a high degree of flexibility, the code, by means 
of input parameters, allows the transform to be varied depending on the region the par- 
ticle is in, the direction it is heading in, and whether the particle is on its first flight. 

A pseudototal macroscopic cross section is defined for each region as follows: 

S *(E) = S t (E)(1 - XK(K) • EXT(K)) 

* 

where E T (E) and S T (E) are the pseudototal and actual total cross sections, respec- 
tively, XK(K) is a direction dependent option which is input for each of the K regions, 
and EXT(K) is a parameter (usually between 0 and 1) which controls the intensity of the 
transform within a given region. 

The quantity XK(K) controls how the transform shall be applied depending on parti- 
cle direction within a region. The allowable options for XK within a region are 

XK = 1, a, y, -a, -/3, -y 

where a, 0, y are the particle direction cosines. If, for example, XK(K) equals 1, this 
signifies that the transform will be applied to all particles within that region independently 
of direction, with an intensity controlled by EXT(K). If XK is chosen to be y, for 
example, then all particles headed in the positive z direction will be biased. A parti- 
cle headed parallel to the z axis (y = 1) will be biased more heavily than one whose 
direction cosine with the z axis is smaller. The value of XK chosen for each of these 
particles will, in fact, be its z axis direction cosine. Particles headed backwards 
(y < 0) undergo no exponential transformation of their flight path. 

The quantity EXT(K) controls the intensity of the transform within a given region, 
independently of particle direction. Setting EXT(K) = 0 is equivalent to not using the 
transform since for this case 


S* = S T (1 - 0) = s T 

* 

Whenever the transform is applied, the pseudototal cross section Srp is used in 

place of the actual cross section 2™ to compute the path length to the next collision. 

1 * 

When EXT is chosen between 0 and 1, Srp is smaller than Erp, and this serves to 
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stretch the path length between collisions. This is the usual case. It should be pointed 
out that EXT can also be chosen to be negative (EXT < 0), in which case S T is larger 
than Erp. This serves to compress the path length, if desired. 

An additional input parameter KBIAS, also allows some control of the exponential 
transform. Setting KBIAS = 0 causes the transform to be applied to all flight paths of 
the particle, while setting KBIAS = 1 causes the transform to be applied only on the 
first flight of the particle, from its birth as a source particle to its first collision point. 

Since the exponential transform is a biasing technique, the bias introduced by its 
use must be removed by suitable modification of the particle weight. If the transform is 
applied within a region and the particle crosses the region without a collision, the bias 
is removed at the boundary crossing by modifying the weight as follows: 


-S 

W' = W e 


[s t (e)-s*(e)] 


If the particle crosses several regions, the weight correction formula is applied at each 
boundary crossing. When the particle reaches the scattering center at a distance D 
within a region, the bias is removed by modifying the particle weight as follows: 

s t (e) -d[s t (e)-s*(e)] 

W’ = W — i e 1 1 

s£(e) 

Histor y termination . - Particle histories are terminated whenever one of the follow- 
ing events occurs: 

(1) The particle's energy is reduced below a specified minimum energy, EMIN. 

(2) The number of collisions in the particle's history exceeds a maximum allowable, 
MAXCOL. 

(3) The particle’s weight is reduced below a specified minimum weight, WMIN. 

(4) The particle escapes into an unbounded region. 

The quantities EMIN, MAXCOL, and WMIN are problem input parameters. When- 
ever the particle weight drops below WMIN, Russian Roulette is played to determine 
whether the history is to be terminated or continued. The Russian Roulette input para- 
meter RR is compared with a random number R. If RR ^ R, the history is continued 
with a new particle weight W’ = W/RR. If RR < R, the history is terminated. 
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Particle Scoring 


Uncollided flux at a p oint detector. - To compute the uncollided flux at a point detec- 
tor, the assumption is made that each particle can be regarded as a microscopic point 
source, emitting l/4n neutrons or photons per unit solid angle. The geometry routine 
computes the total number of mean free paths along the flight path and the uncollided flux 
from each particle becomes 


_ We' A 
^unc 9 
4jkT 

where W is the weight of the particle, X is the total number of mean free paths between 
particle and detector, and d is the distance between particle and detector. Because of 
these assumptions, the uncollided flux calculation has no meaning for a monodirectional 
or anisotropic source. 

Scattered flux at a point detector . - In addition to the uncollided flux, the code es - 
timates the scattered flux at each detector point by computing the probability at every 
collision that a particle will scatter in a direction so as to be headed towards the detec- 
tor and will arrive there with no further collisions. 

The equation used to estimate the scattered flux at a point detector from each col- 
lision is given by 


— tt-.E) 


^sc 


W dfi 
2 


a(E) 


where W is the particle weight after collision, d is the distance from collision point 
to detector, dcr/dft^, E) and cr(E) are, respectively, the differential and total scattering 
cross sections for a particle of energy E, and X is the total number of mean free paths 
from the collision point to the detector. The ratio (da/dS2)/a represents the probability 
of a particle scattering through an angle \f/ in the laboratory system, so as to be headed 

o 

towards the detector; 1/d is the solid angle subtended by a unit sphere at the detector; 
and finally, e -Al is the probability of the particle arriving at the detector with no further 
collisions. The flux estimator is computed and stored for each detector as a function of 
the particle energy and also, if desired, the polar and azimuthal angles of arrival at 
the detector. 

The computation of the scattered flux at a point detector begins with the calculation 
of the polar scattering angle in the laboratory system. Let the spatial coordinates of a 
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point detector and a particle be x., y., z., and X, Y, and Z, respectively. 

X liu X 

distance between the particle and the 1 1 detector is given by 


d A = jpC - x.) 2 + (Y - y.) 2 + (Z - z ;.) J 


11/2 


and the direction cosines of the flight path to the detector are 


“i = 




y i = 


z. 

1 


- Z 


If the direction cosines of the particle before collision are a, /3, and y, then the labora- 
tory system polar scattering angle through which the particle must scatter in order to be 
headed toward the detector is given by 


cos i// = ao'j + /3/3j + yy ^ 

Turning first to the treatment for gamma rays, the angular scattering probability 
in the laboratory system is next computed from the Klein-Nishina formula by use of the 
expression 


— (xfy E) 2 <3 9 9 

dS2 = r o (P + - P* sim £) 

ct c (E) 2 ct c (E) 


where 


P = 


1 


1 + — — — (1 - cos \f/) 
0. 511 
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r is the classical electron radius, and a (E) is the gamma-ray scattering cross sec- 
o c 

tion per electron at energy E. The cross section cr„(E) is computed from the formula 


<7 c (E) = 2,1-2 



- 2 Q (1 t . Q) - In (1 + 2Q) +J-ln(l + 2Q) 
1 + 2Q _ 2Q 


1 + 3Q 
(1 + 2Q) 



where Q = E/0. 511 for the gamma -ray energy in MeV. These equations are computed 
using the gamma -ray energy after scattering which is given by 


1 + — - — (1 - cos xjy) 
0. 511 


where E is the gamma -ray energy before collision. 

The treatment is not as simple for neutrons as for the gamma rays. A conversion 
must be made from the laboratory scattering angle ;// to the center -of -mass scattering 
angle X to find the angular scattering probability. The angular scattering probability, 
computed in the center -of -mass system, must then be reconverted to the laboratory 
system. The center -of -mass scattering angle X is given by 


cos X 


cos 2 l// 


,Sc 


1 -f cos cost + A^B - 1 
_ 41*1/3 


where A and B have been previously defined. (See Collision Processes section. ) 

The angular scattering probability in the center -of -mass system is next computed. 
For those cases in which elastic scattering is isotropic in the center -of -mass system 
and for all inelastic scattering events, the probability of scattering per unit solid angle 
into the angle X is l/4n. 

When neutron elastic scattering is anisotropic in the center -of -mass system, the 
probability of scattering per unit solid angle into the angle X is computed from the 
cumulative probability distribution tables (library 10 data) for the element of the colli- 
sion. 

The scattering probability in the center -of -mass system is then converted to the 
laboratory system by means of the equation 

— O^E) — (X,E), 2 \3/2 

dft dfi (,A + 2A cos A + 1 ) 

a c (E) a c (E) a 2 |A + cosa| 
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The equation used to compute the neutron energy after scattering is the same as that 
used in the collision processes, namely, 

E' = - (1 + A 2 B + 2AB 1//2 cos A) 

(A + 1) 2 


As in the collision processes, hydrogen scattering allows for a much simpler calcu- 
lation. For hydrogen, the scattering probability per unit solid angle in the center -of - 
mass system is isotropic. Therefore, 


— <*, E) 
dfi 

a c (E) 


I2 3 / 2 (1 + cosA) 1 / 2 

47T 


and 


cos A = 2 


cos 2 ip 


Combining these, we obtain for hydrogen 


da 


(*,E) 


dfi 

_ a c (E) 


cos ip 

n 


o 

with the special case E' = E cos for hydrogen. When scattering in hydrogen is in- 
volved, if the laboratory scattering angle \ p is larger than 90° (cos xp < 0), the calcu- 
lation of the flux at the detector is terminated since neutron scattering in hydrogen 
cannot physically occur through angles larger than 90°. 

Region averaged fluxes . - The flux per unit source strength for a region is obtained 
by estimating the expected track length per unit source strength through the region 
divided by the volume of the region. The track lengths in each region are summed for 
all particles that can reach the region. The expected track length per unit source parti- 
cle is the sum of the track lengths divided by the number of source particles. 

The expected track length through a region J for a particle is given by 


TLj = W 



p(x)x dx + W 


S e 


-SS 
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where TLj is the expected track length through region J, W is the weight of the par- 
ticle, S is the distance across region J, S is the total macroscopic cross section for 
region J, and p(x) is the probability per unit distance of a particle making a collision 
at x. 

The first term of the preceding equation gives the expected track length contribution 

for those particles traveling a distance x less than S, and the second term gives the 

expected track length contribution for region J for those particles that travel beyond 

- y)x 

it. Since p(x) = 2 e , the equation becomes 


TL J = W 


/.S 

I Sx e _Sx dx + WS e~ SS 

J o 


w 

s 


1 - e" SS (1 + SS)1 +WS e _SS 


W 

s 


1 - e _SS (1 + SS - 2S)1 = — (1 - e _SS ) 

J 2 


When the material in the region is a vacuum the expected track length through the 
region is simply 


TL 


vac 


= W • S 


The weight of the particle entering a region not containing the collision point is re- 
defined at the boundary as 


W’ = W e -EXP 

where EXP is the sum of the mean free paths traveled by the particle since the last 
boundary. 

Particle leakage tape and p rint out table . - Particles leaking out of the system may 
be tracked in a subsequent problem, if desired, by generating a leakage tape. Also, a 
leakage printout table that sorts the particles leaking out according to energy and angle 
of emission with the Z-axis'can be obtained. These options are available only if region 
averaged fluxes are requested. (For details, see the section DATA INPUT.) 

The spatial location, direction cosines, energy, weight, and region number of the 
particle at the time of its escape are recorded on the leakage tape. Since these are the 
same parameters written by the primary particle source generator, this tape can be 
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used in a subsequent run and the leakage particles can be tracked. The leakage printout 
table can sometimes be used to obtain the angular and energy distribution for a point 
source representation of the system. 

Energy deposition . - Energy deposition for neutrons and photons is computed at each 
collision whenever the incident particle energy E is greater than an input parameter 
(ELD) . The energy deposition from a neutron or photon absorption is computed as 

ED(1) = W E(1 - P (E)) 
s 

where W is the particle weight before collision, E is the particle incident energy, 

P (E) is the scattering probability for particles of energy E, and 1 - P (E) is the ab- 
s s 

sorption probability for particles of energy E. 

The energy deposition from neutron or photon scattering is computed by use of the 
formula 


ED(2) - W P (E)(E - E’ - EJ 

S X 

In order to use the same formula for energy deposition from scattering for both 
neutrons and photons, the quantities E’ and E„ must be defined differently depending 
on the type of particle and type of scatter event. 

For neutrons, E’ is defined as the neutron energy after scattering. For neutron 
elastic scattering E x = 0, but for neutron inelastic scattering E x is the excitation 
energy of the nuclide. The excitation energy E v is given off in the form of inelastic 
gamma rays, and the energy deposition from these should be treated as a separate pro- 
blem. 

For gammas, if the event is a Compton scatter, E x is zero, and E' is the photon 
energy after scatter. If the event is a pair production, then E’ = E x = 0. 511 MeV, 
since, for a pair production event, the original photon disappears and is replaced by two 
photons, each of energy 0. 511 MeV. 

Statistical analysis . - The program calculates the standard deviation a and the 
percent deviation for (1) energy deposited in regions, (2) fluxes at point detectors, and 
(3) track length fluxes in regions. The method used in the code is that of reference 6. 
The random variables (energy deposition, fluxes, track lengths) are tallied as group 
scores accumulated from tracking a group of n histories. The group size n is an in- 
put parameter chosen by the user for the N histories in which G = N/n groups are 
processed in the problem. The value of n should be chosen so that G > 10 for best 
statistical results. 

The basic equation used in the calculation of the standard deviation per history is 
given by 
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where 


<r-l. 1 


n \l gT 


V s g -( s ) : 


G 

2 >, 


g 


_ z 

s 2 = §C-L- 

* G 


g 


The percent deviation (PC DEV) of each random variable x is the standard devia- 
tion divided by the average. 


PC DEV = 100 = 
X 


where X is the calculated average value of x. 


History Tracking Techniques 

During particle tracking, the code makes four passes through the geometry to de- 
termine (1) distance to the next collision, (2) estimated track lengths in regions, and 
(3) uncollided and (4) scattered fluxes at point detectors. Since the program repeats 
(1), (2), and (4) at each collision, a significant amount of computing time is spent in 
geometry calculations. The program has built-in features to aid in reducing geometry 
tracking time. 

As the code calculates the distance to the next collision, track length estimation and 
particle tracking proceed together up to the collision point. The track lengths in the re- 
maining regions beyond the collision point in the direction of travel are then calculated 
before the code resumes history tracking. 

Geometry tracking can be reduced by terminating the track length estimation when 
the contribution from that collision becomes negligible. The code provides a cutoff 
option to accomplish this. When the number of mean free paths of the track length is 
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greater than the input cutoff parameter (TLCUT), the track length calculation is ter- 
minated and history tracking resumed. 

A cutoff may also be used to terminate the estimation of fluxes at point detectors 
when these become negligible. The input parameter EPC(K)is a detector dependent quan- 
tity, which may vary for each of the K detector points. When the cumulative number of 
mean free paths from the collision point to the K m detector is greater than EPC(K), flux 
estimation to this detector is terminated, and the code proceeds to the next detector. 

The setting of these cutoff values are at the complete control of the user. 


SECONDARY GAMMA SOURCE GENERATOR 

The function of the secondary gamma source generator is to generate secondary 
gamma photons resulting from inelastic scatter (n, n',y) or from neutron capture (n, y) 
reactions with the shield materials and to write a secondary gamma source tape. The 
main features are summarized as follows: (1) Secondary gammas are generated for 
only two interactions, inelastic scatter n, n’,y or neutron captures n, y. (2) The 
routine randomly selects only one gamma ray from a capture gamma spectrum, but cal- 
culates a weight that preserves the total energy released per capture. (3) The code 
analyzes collisions that give rise to inelastic or capture gammas in a single computer 
run. The code requires, as input for the secondary gamma -ray analysis, a collision 
parameter tape and a card input deck for the neutron capture gamma spectrum for each 
element of interest. Card input is not required to generate gamma rays resulting from 
inelastic scatter. 

The output is a gamma -ray source tape to be used as input for a subsequent analysis 
of the secondary gamma sources. Written on the tape are the spatial coordinates x, y, 
and z, direction cosines a, 13, and y, energy E, weight W, and region number NR of 
each secondary gamma -ray source particle. A code option INEL, allows the user to get 
inelastic and capture gamma sources on separate tapes, or both types of gamma -ray 
sources on a single tape. 

Generation of secondary gammas from inelastic scatter is the simpler of the two 
methods. The initial energy E and final energy, E’ of the inelastically scattered neu- 
tron are taken from the collision tape. The gamma -ray energy E^ is simply the dif- 
ference in these energy values: 


E = E - E f (1) 

r 

The weight of the gamma ray W^ is equal to the weight of the incident neutron times the 
probability that no absorption took place. Thus 
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I 



( 2 ) 


W = W P (E) 
y ns'' 

where P (E) is the scattering probability at energy E. The spatial coordinates x, y, 
s 

and z and the region of birth of the gamma ray (NR) are taken directly from the colli- 
sion tape. 

The angular distribution of the secondary gamma rays is assumed to be isotropic in 
the laboratory system. Consequently the initial direction cosines a, p, and y are 
chosen at random by use of the rejection technique illustrated in figure 2. 

The geometric parameters for capture gamma sources (x, y, z, a,/3,y, NR) are ob- 
tained the same way as outlined for inelastic gamma sources. The calculations of the 
gamma -ray energy E^ and weight however, are more involved. 

The chief input data for generating capture gammas is a table which gives the num- 
ber of photons Nj j at gamma energy E^, emitted when a neutron of incident energy 
E^ is captured by element j. When a neutron of energy E^ is captured in element j 
the total energy released Q(Ej) is 

Q(E i> = E N J, A W 

k 


For a given element, the probability ^ 
ISk is 


of picking a single gamma photon of energy 


l, k 


A A 3, A 

Q(Ej) 


IX, A 


(4) 


The table N. , , is input for each element that can produce capture gammas. 

Then Q(Ep is computed internally for each neutron incident energy and element. The 
probability table P , v is then constructed for each element and further modified to a 

4 3 K 

cumulative probability distribution k . The element in which the absorption occurs 
is available from the neutron collision tape. The energy of the secondary gamma photon 
E k is chosen from the cumulative probability distribution by obtaining a random number 
Rj and choosing the smallest value of k such that 

IX A *1 (5) 

k 
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The number of gamma rays n* of energy required to preserve the total energy per 
capture Q(Ep is 


Q(E,) 

n’ = — 

E k 


( 6 ) 


Only one photon of energy E^ is created, but the factor n’ is used in determining the 
weight of the photon. 

The weight of the photon created is calculated as the product of the neutron 

weight W , the probability of an absorption l"l - P (E)~|, where P (E) is the scattering 
n s s 

probability and the total energy per capture weighting factor n’. Thus 

W y - W „& - P S <E)>’ (V 

The absorption probabilities are calculated from the microscopic cross sections used 
in the history tracking analysis. 

This procedure generates capture gammas at every collision point where the ab- 
sorption probability is nonzero. This can result in the generation of an excessive num- 
ber (from a computer time viewpoint) of secondary gammas originating from a relatively 
few primary neutrons. 

At the user's option capture gammas can be generated either at every collision 
where the absorption probability is nonzero or only for those collisions where the ab- 
sorption probability jjL - P g (E)] is greater than a random number. For this second op- 
tion, the absorption probability is compared with a random number at every collision 
point. If the probability is less than the random number, no capture gamma is generated. 
If the probability is greater, a gamma is generated. Its energy is obtained as described 
previously, but its weight is given by 

W = W n' 
y n 

This procedure, which can be regarded as a form of Russian Roulette, yields a 
smaller number of secondary gammas and is more economical of computer time. 


PREPARATION OF INPUT DATA 
CODE DESCRIPTION 


The code is written completely in FORTRAN -IV for the IBM 7094 computer. The 
overlay feature has been used to allow sufficient storage for large problems; variable 


33 



dimensioning has been incorporated. The program consists of a main link and five de- 
pendent links (see fig. 1). At any given time, the main link and only one of the depend- 
ent li nk s are stored in core; the remaining links are stored on the overlay tape. 

The main link acts as an executive routine, calling in the other links as required to 
manipulate data or perform calculations. The first link reads in and prints out input 
data and cross sections and computes the storage requirements of the problem. The 
second link calculates primary particle source parameters from spatial, angular, and 
energy distributions read in by the user. These parameters are written on a tape and 
used by the particle tracking routines. If secondary gammas or leakage particles are 
being tracked from a previous problem, then a physical source tape will be available 
for this run. In this case the second link is bypassed with control proceeding directly 
from the first link to the third link. 

The great majority of computing time is spent in the third link which actually tracks 
and scores the particles. After the desired number of particles has been followed, con- 
trol passes to the fourth link which prints out the problem results. 

If no secondary gamma source is to be written, the problem terminates at this 
point; otherwise control proceeds to the fifth link, which generates a secondary gamma 
source tape. When this link is completed the problem ends. 

The last random number is punched out at the end of each run and is available for 
use in starting a new sequence of random numbers for subsequent problems. 


DATA INPUT 

The code requires that all input data defining the problem be read in on cards at 
execution time. A summary of the read order for the data follows: 

(1) Read in the problem control parameters. 

(2) Read in the additional problem input consisting of atomic weights, atom densities, 
printout bounds for energy, polar and azimuthal angles, and point detector data. 

(3) Read in the input library data. These include libraries 1, 2, and 7 which give the 
boundary descriptions, region descriptions, and exponential transform instructions, 
respectively. Also, required is one library 6 for each element, and for neutron pro- 
blems libraries 9 and 10 are required for those elements where inelastic scatter or 
anisotropic elastic scatter occurs. For the user’s convenience, the library data may 
be read in any order. For example, library 6 for a given element can be read first, 
followed by libraries 1, 7, etc. , etc. 

(4) Read in the primary source generator control parameters. 

(5) Read in remainder of primary source generator input. 
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Items (4) and (5) are not required if an external source tape from a previous run 
is supplied. 

(6) Read in response function data, where desired. 

(7) Read in capture gamma input data, for those problems where a secondary gamma 
source tape is to be generated. 

The format for preparing a problem data deck is given in tables I to XII. 

Problem control parameters . - The problem control parameters, consisting of 
fixed and floating point data, are required input for every problem. The fixed point 
parameters are used in determining the amount of storage required for the problem. 

The floating point and fixed point parameters are both used as limiting controls on vari- 
ous aspects of the problem during execution. The format for preparing these data are 
given in table I. 

Problem description data . - The atomic weights and atomic densities are required 
to describe the material of the problem. Also, the energy and angular bin values for 
the printout of the energy and angular fluxes must be read in. Table n gives the format 
for preparing these data. 

Library 1 . - Library 1 contains information describing boundaries for the problem 
geometry. The formats to be used are given in table in. Library 1 is necessary for 
every problem. 

Library 2 . - In library 2 the various regions in the geometry are described by 
giving the number of boundaries and the sign of each boundary that defines the re- 
gion, along with the material number and the volume of the region. One library 2 is 
necessary for every problem. Table IV shows the format to be used in preparing a li- 
brary 2. The material within any region must be the same for all points within the re- 
gion. 

Usually, the first region described is the source region. Regions may be listed in 
any order desired, but other arrangements may prove to be more economical. For com- 
puter economy the number of regions in the geometry should be kept to a minimum, con- 
sistent with accuracy of representation and type of output desired. 

Lib rary 6 . - A library 6 is required for every element in the problem. This library, 
for which the input format is given in table V, lists, for arbitrary energy values, the 
total, scattering (elastic plus in-elastic), and elastic scattering cross sections for neu- 
trons or the total, Compton plus pair production, and Compton scattering cross sections 
for gamma rays. The energy values must be the same for every element in the problem. 
The cross sections are listed for the highest energy first and then in descending order to 
the lowest energy. The library 6's may be loaded behind the problem input data in any 
order convenient to the user. 
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TABLE I. - PROBLEM CONTROL PARAMETERS 


Card 

Format 

Card 

Input 

Description 

group 


column 

item 


1 

12A6 

1-72 

Title 

Any alphanumeric identification in card 





columns 1 to 72 

2 

018 

1 - 18 

RMN 

Starting random number 

3 

1415 

1 - 5 

IDP 

Problem ID number (any fixed point num- 





ber for user’s convenience) 



6-10 

NRG 

Particle type: NRG = 0, neutrons; 





NRG = 1, gammas 



11 - 15 

NGPSZ 

Number of particles per group (for 





statistical grouping) 



16 - 20 

NEL 

Number of elements in problem 



21 - 25 

MAT 

Number of materials in problem 



26 - 30 

NREG 

Number of regions in problem 



31 - 35 

NBD 

Total number of boundaries 



36 - 40 

NBDMX 

Maximum number of boundaries in any 





given region 



41 - 45 

INLIBR 

Number of libraries to be read in (The code 





compares this number with number of 
libraries actually supplied and with a 
number internally calculated to be 
necessary. If the three numbers are 
not the same, an error message will 
be printed out. ) 



46 - 50 

NGY 

Number of energy values in library 6 cross 





section tables 



51 - 55 

NXSEC 

Cross-section printout option: NXSEC = 0, 





does not print out cross sections; 

NXSEC = 1, prints out cross-section data 



56 - 60 

ND 

Number of point detectors: ND = 0, does not 





calculate flux at point detectors; ND = K, 
calculates flux at K point detectors 



61 - 65 

NTLF 

Region flux option: NTLF = 0, does not cal- 





culate fluxes in regions; NTLF = 1, calcu- 
lates fluxes in regions (by track length 
method); NTLF = 2, same as NTLF = 1, 
but also prints leakage data; NTLF = 3, 
same as NTLF = 2 but also writes leakage tape 



66 - 70 

LTNO 

Tape number to be assigned to leakage tape (ig- 





nore unless NTLF = 3. ) 



TABLE I. - Continued. PROBLEM CONTROL PARAMETERS 


Card 

Format 

Card 

Input 

Description 

group 


column 

item 


~T~ 

915 

1- 5“ 

IEP 

Number of output energy bins 



6 - 10 

IGP 

Number of output polar angle bins 



11 - 15 

IPP 

Number of output azimuthal angle bins 



16 - 20 

NSOR 

Region number in which source is located 



21 - 25 

ISTAP 

Source tape option: ISTAP = 0, source particles 





are generated internally by primary particle 
source generator; ISTAP = K, source par- 
ticle data are available on tape number K 



26 - 30 

MAXCOL 

Maximum number of collisions per particle 





allowed (Any particle undergoing a larger 
number of collisions will have its history 
terminated. ) 



31 - 35 

KBIAS 

Exponential transform on first flight option: 





KBIAS = 0, exponential transform applied 
(as determined by library 7) for all flight 
paths of the particle; KBIAS = 1, exponen- 
tial transform applied only for first flight 
path of particle (No transform applied for 
succeeding flight paths) 



36 - 40 

ISPON 

Response function option: ISPON = 0, does not 





compute response function results; ISPON = K, 
computes results for (0 < K < 8) response 
functions entered 



41 - 45 

NSEG 

Collision tape write option: NSEG = 0, does not 




i 

write collision parameters on tape; NSEG = 1, 
writes collision parameters on scratch tape for 
use in secondary gamma generator 

5 a 

615 

1 - 5 

NL9 

Number of elements that require library 9 data 



6-10 

INE9 

Maximum number of incident neutron energies 





in any library 9 for this problem 



11 - 15 

IEX 

Maximum number of nucleus excitation levels 





in any library 9 for this problem 



16 - 20 

NL10 

Number of elements that require library 10 data 



21 - 25 

1NE10 

Maximum number of incident neutron energies in 





and library 10 data 



26 - 30 

DPR 

Maximum number of cosines of scattering angles 





in any library 10 for this problem 


a Omit card 5 for gamma problems. 



TABLE I. - Concluded. PROBLEM CONTROL PARAMETERS 


Card 

Format 

Card 


Description 

group 


column 



6 

8E9.4 

1 - 9“ 

HIST 

Number of particle histories to be tracked 



10 - 18 

EMAX 

Maximum energy (MeV) for any particle or 





cross-section table in this problem 



19 - 27 

EMIN 

Minimum energy (MeV) for this problem; also. 





cutoff energy for any particle, below which 
its history will be terminated 



28 - 36 

WMIN 

Weight cutoff for particle 



37 - 45 

ELD 

Energy deposition cutoff value (MeV); energy 





deposited below this value will be assumed 
to be negligible 



46 - 54 

ETM 

Thermal energy value (MeV); neutrons below 





this value scatter isotropically without 
energy loss (Ignore for gammas. ) 



55 - 63 

RR 

Russian Roulette parameter (Russian Roulette 





played with particles whose weight is less 
than WMIN; RR is their chance of surviving. ) 



64 - 72 

EPSL 

Small incremental distance used to move a par- 





tide’s position off the boundary in geometry 
calculations 

7 

1 2E9.4 

1-9 

CTIME 

Cutoff time parameter (The clock is read at the 





completion of every batch of 100 histories, and 
the time remaining is computed (Estimated 
time for problem minus elapsed time.) If the 
time remaining is greater than CTIME, another 
batch of 100 histories will be tracked. If not, 
tracking is abandoned and the output is compiled 
for those histories already completed. CTIME 
is a device for obtaining at least partial results 
when a problem runs longer than the time re- 
quired to track 100 histories. Except for this, 
no specific recommendations are given. Try 
CTIME =1.0 min for a start. ) 



10 - 18 

TLCUT 

Track length cutoff parameter (If a particle’s 





cumulative number of mean free paths through 
the geometry exceeds TLCUT, the track length 
calculation is abandoned for that particle. ) 



TABLE II. - PROBLEM DESCRIPTION DATA 


Card 

group 

Format 

Card 

column 

Input 

item 

Description 

1 

(8E9.4) 

1 -72 a 

ATWT(I) 
(1= 1,NEL) 

Atomic weights of elements (For 
neutrons, the atomic weights of 
elements requiring libraries 9 
and 10 should be read in first. ) 

2 

(8E9. 4) 

1-72 

DN(J, I) 
(J= 1, MAT) 
(I = 1, NEL) 

Atomic density of element I in ma- 
terial J. All the densities for 
elements of material 1 are read 
in first, followed by all the densities 
for the elements of material 2, etc. 

3 

(8E9.4) 

1-72 

ESI(J) 

(J= 1, NEL) 

Energy below which elastic scattering 
of neutrons with element J is as- 
sumed to be isotropic in the center— 
of -mass system (Omit for gammas.) 

4 

(8E9.4) 

1-72 

i 

EPRINT(K) 
(K = 1, IEP+1) 

Energies bounding the output energy 
bins (Read-in MeV and, in ascending 
order, the first energy bound should 
be EMIN and the remaining values 
the upper bounds of each of the output 
energy bins. ) 

5 

(8E9.4) 

1-72 

GPRINT(L) 
(L = 1, IGP+1) 

i 

Angles bounding polar angle output bins 
(Read-in degrees and in ascending order. ) 

b 6 

(8E9.4) 

. .. 

1 1-72 

PPRINT(L) 
(L = 1, IPP+1) 

Angles bounding the azimuthal angle output 
i bins (Read in degrees and in ascending 
■ order. ) 

C 7 

(8E9.4) 

1-9 

37-45 

10-18 
46-54 
19-27 
55-63 
28 -36 
64-72 

EPC(K) 

XD1(K) 

YD1(K) 

ZD1(K) 

Cutoff parameter for detector K, in mean 
free paths 

Rectangular spatial coordinates for 
detector K 


a l -72 implies that data are read according to the format in 8 fields, 9 columns 


wide, for the entire card. 
b Omit cards 6 and 7 if ND = 0. 

c Read order: Two sets of detector parameters per card as follows: 

EPC(K), XD1(K), YD1(K), ZD1(K), EPC(K + 1), XD1(K + 1) YD1(K + 1), ZD1(K + 1). 



TABLE IB. - LIBRARY 1 INPUT DATA FORMAT 


Card 

Format 

Card 

Input 

Description 

group 


column 

item 


1 

219 

1 - 9 

LIBT 

Library type (LIBT = 1 for library 1) 



10 - 18 

I 

Number of boundaries describing 
the geometry (I = NBD) 

a 2 

(19, 7E9. 4) 

1-9 

IBT(K) 

Boundary type; IBT(K) = 1, general 





quadratic equation; 1BT(K) = 2, 
spheres, cylinders; IBT(K) = 3, 
planes perpendicular to the X-axis; 
IBT(K) = 4, planes perpendicular 
to the Y-axis; IBT(K) = 5, planes 
perpendicular to the Z-axis 



10 - 18 

AF(K) 

Coefficients for boundary K 



19 - 27 

BF(K) 




28 - 36 

CF(K) 




37 - 45 

XF(K) 




46 - 54 

YF(K) 




55 - 63 

ZF(K) 




64 - 72 

DF(K) 


b 2(a) 

(19, 3E9. 4) 

1 - 9 

IBT(K) 

Boundary type coefficients for boundary K 



10 - 18 

EF(K) 




19 - 27 

FF(K) 




28 - 36 

GF(K) 



a One card 2 is required for each boundary. 

^Card 2(a) is needed only when a type 1 boundary is specified. It immediately 


follows card 2 for the same boundary. 



TABLE IV. - LIBRARY 2 INPUT DATA FORMAT 


Card 

group 

Format 

Card 

column 

Input 

item 

Description 

1 

15 

1-9 
10 - 18 

LIBT 

I 

Library type (LIBT = 2 for library 2) 
Number of regions in problem 
(I = NREG) 

a 2 

* 

1-5 

NB(J) 

Number of boundaries defining region J. 
(Max. of 10). Negative sign on NB(J) 
indicates that region J is an outside 
region 


14 

6-9 

MATREG(J) 

Material number in region J 


E9.4 

10 - 18 

V(J) 

Volume of region J in cm^ (Enter any 
nonzero number unless region fluxes 
are being computed for region J, in 
which case the correct volume is 
required. ) 


6(15, 14) 

Begin 
19 - 23 ' 

Begin 
24 - 27 

IBN(J,K) 
MPR(J, K) 

K^ 1 boundary surrounding region J. 

(K = 1, NB(J)). (Note: Boundary 
numbers must be given a plus or 
minus sign that corresponds to the 
direction index as discussed in sec- 
tion Geometry Method. ) 

Most probable region of entry across 
the K th boundary from region J 


a One card 2 is required for each region. The reading order is as follows: 
(NB(J), MATREG(J), V(J), IBN(J, 1), MPR(J, 1), IBN(J, 2) MPR(J, 2). . . 
IBN(J, K), MPR(J,K)). 


TABLE V. - LIBRARY 6 INPUT DATA FORMAT 


Card 

group 

Format 

Card 

column 

1 

319 

1 - 9 
10 - 18 



19 - 27 

2 

8E9.4 

1 - 9 
10 - 18 
19 - 27 
55 - 63 



28 - 36 
64 - 72 


Input 

item 

LIBT 

I 


l ENERGY(L) 
TM(J, L) 
SM(J, L) 


ESM(J,L) 


Description 


Library type (LIBT = 6) 

Number of energy values at which 
cross sections are read in (I = NGY) 
Element number (elements are* numbered 
in same order that their atomic weights 
are listed) 

Energy value in MeV 

Total microscopic cross section for element J 
Microscopic scattering cross section for ele- 
ment J (elastic plus inelastic scattering 
cross section for neutrons or Compton 
plus pair- production cross section for 
gamma rays) 

Microscopic elastic scattering cross section 
for neutrons or Compton scattering cross 
section for gammas for element J. 


a Reading Order: 
TM(J, L + 1), 


(ENERGY(L), TM(J,L), SM(J,L), 
SM(J ? L + 1), ESM(J,L + 1)). 


ESM(J,L), ENERGY(L + 1), 



The units used for library 6 cross sections must be consistent with the units speci- 
fied earlier for the atomic density DN(J,I) of element I in material J, since the product 
of a library 6 cross section and DN must have the units of inverse length. Two con- 
sistent sets that satisfy this requirement occur if the library 6 cross section is in square 
centimeters per gram and DN is in grams per cubic centimeter, or if the library 6 
cross section is in barns per atom and DN is in atoms per barn- centimeter. 

Library 7 . - One library 7 is necessary for each problem. Library 7 describes 
how the exponential transformation will be applied in each region of the problem and the 
biasing parameter used to obtain a pseudo cross section which is used in the sampling of 
path lengths. The input format for preparing a library 7 is shown in table VI. 

Library 9 . - This library is omitted for gamma-ray problems. For neutron pro- 
blems, a library 9 is required for every element where inelastic scattering can occur. 
The library lists the cumulative probability of exciting each energy level for various in- 
cident neutron energies. The format for preparing a library 9 is shown in table VII. 

Library 10 . - A library 10 is required for each element in which neutron elastic 
scattering is anisotropic in the center -of -mass system. If the energy values ESI(J) 
for the element J is less than the maximum energy EMAX, a library 10 is required for 
this element. A library 10 is not required for neutron elastic scattering by hydrogen 
or for gamma -ray problems. 

Library 10 lists the cosines of the center -of -mass scatter ing angles that divide the 
integral of the differential elastic scattering cross section over the solid angle into equal 
parts. The cosines are listed for various incident neutron energies. The format for 
preparing library 10 data is shown in table VIH. 

Primary source control parameters . - The controls are required to generate the 
source particle parameters for the option specified by the user. The format for pre- 
paring these input data is shown in table IX. 

Primary source generator input . - The source input requires a set of relative power 
distributions to completely describe the spatial (rectangular, cylindrical, or spherical) 
distributions of the source (neutrons or gammas). Distributions giving the fraction of 
particles in each energy or angular interval are also required to define the energy and 
angular distributions of the primary source. The format for preparing these data are 
given in table X. 

Response function data . - Table XI gives the input format for the response function 
data. This table is omitted when ISPON = 0. For ISPON = N (0 < N < 8), N sets of 
cards are read in. The first card identifies the response function, and the rest the 
actual data. It should be emphasized that the response function data are entered as 
averages over the output energy bins. 
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TABLE VI. - LIBRARY 7 INPUT DATA FORMAT 

Input Definition 

item 

LIBT Library type (LIBT = 7) 

I Number of regions (I = NREG) 

NEXT(M) Exponential transform direction op- 
tion for region M: 

NEXT(M) = 1, XK = 0 (transform not used) 
NEXT(M) = 2, XK = 1 (transform applied 
independently of direction) 

NEXT(M) = 3, XK = a 1 (transform applied to 
NEXT(M) = 4, XK = -aj x or -x component) 
NEXT(M) = 5, XK = /3 1 (transform applied to 
NEXT(M) = 6, XK = -/3J y or -y component) 
NEXT(M) = 7, XK = y 1 (transform applied to 
NEXT(M) = 8, XK = -yj z or -z component) 
EXT(M) Exponential transform intensity parameter 
for region M (0 < EXT(M) < 1. 0) 


a Reading Order for card 2 is as follows: (NEXT(l), EXT(l), NEXT(2), EXT(2). . . . 
NEXT(NREG), EXT(NREG). 

^XK and EXT are used to calculate a pseudo cross section by the following 
* 

expression: £(E) = E„(E)(1 - XKT(M)). 

T 1 
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TABLE VII. - LIBRARY 9 INPUT DATA FORMAT 


Card 

group 

Format 

Card 

column 

Input 

item 

Description 

1 

419 

1-9 
10 - 18 

19 - 27 

28 - 36 

LIBT 

NMESH(J) 

J 

NEXP(J) 

Library type (LIBT = 9) 

Number of incident neutron energies 
entered for element J 
Element number (Elements are numbered 
in same order as atomic weights. ) 
Number of nuclear levels entered for 
element J 

2 

8E9.4 

1-72 

EEK(J, N) 

(N = 1, NEXP(J)) 

Excitation levels for element J in MeV 
(listed in ascending order, first 
excitation level should be zero) 

3 

8E9.4 

1-72 

EMESH(J,K) 

(K = 1, NMESH(J)) 

Incident neutron energies for element 
J in MeV (These energies may differ 
for each element and are listed in 
decending order. ) 

4 

8E9.4 

1 - 72 

SEP(J, K, N) 

Cumulative probability of incident 
neutron of energy K exciting energy 
level N of element J (Cumulative 
probabilities are entered for every 
excitation level for the first neutron 
incident energy, then for every ex- 
citation level for the second incident 
energy and so on. ) 



TABLE VIA. - LIBRARY 10 INPUT DATA FORMAT 


Card 

Format 

Card 

Input 

Description 

group 


column 

item 


1 

419 

1 - 9 

LIBT 

Library type (LIBT =10) 



10 - 18 

IPROB (J) 

Number of cosines of the scattering 





angle listed for element J 



19 - 27 

J 

Element number (elements are num- 





bered in same order as atomic weights. ) 



28 - 36 

NCROS(J) 

Number of neutron incident energies entered 





for element J 

2 

(8E9.4) 

1-72 

ESCAT(J, N) 

Incident neutron energies (MeV) (listed in 




(N = 1, NCROS(J)) 

descending order from 1 to NCROS(J) 
for element J 

3 

(8E9.4) 

1-72 

DES(J,N,K) 

Cosines of the center -of -mass scattering 




(N = 1, NCROS(J)) 

angle for a neutron of incident energy N 




(K = 1, IPROB(J)) 

with element J; IPROB cosine values are 
entered, in descending order from +1 
to -1, for the first incident neutron energy, 
then another set of IPROB cosine values 
for the second incident neutron energy, 
and so on 



TABLE IX. - PRIMARY SOURCE CONTROL PARAMETERS 


Card 


Format 


Card 


group 


1 


1015 


column 

1-5 


6-10 


11 - 15 


16 - 20 


IC 


ID 


IE 


NX 


Description 


Source geometry option: 

IC = 0, point source 

IC = 1, rectangular parallelepiped source 

IC = -1, rectangular parallelepiped source with 

spatial biasing 

IC = 2, cylindrical source 

IC = -2, cylindrical source with spatial biasing 
IC = 3, spherical source 

IC = -3, spherical source with spatial biasing 
Source angular distribution option: 

ID = 1, monodirectional source 
ID = 2, isotropic angular distribution 
ID - -2, isotropic angular distribution with polar 
and azimuthal angular biasing 
ED = 3, anisotropic angular distribution 
Source energy option: 

IE = 1, monoenergetic source 
IE - 2, energy spectrum 

IE ss -2, energy spectrum with energy biasing 
Number of spatial intervals of source in first coordinate 
direction (Ignore for IC = 0); see table below to de- 
termine which is first, second, etc. , direction 
depending upon source geometry 


Source 

Coordinate direction 

geometry 

First 

Second 

Third 

Rectangular 

X 

y 

z 

Cylindrical 

r 

e 

z 

Spherical 

r 

0 

0 


21 - 25 
26 - 30 
31 - 35 
36 - 40 


NY 

NZ 

NSGY 

NANP 


41 - 45 


NANZ 


Number of spatial intervals of source in second coordinate 
direction (Ignore for IC - 0) 

Number of spatial intervals of source in third coordinate 
direction (Ignore for IC = 0) 

Number of energy intervals defining source spectrum 
(Ignore for IE = 1) 

Number of polar angular intervals used to describe 
the source angular distribution (used only for 
ID = -2 or ID = 3) 

Number of azimuthal angular intervals used to describe 
source angular distribution (Used only for ID = -2 or 
ID = 3) 


46 - 50 


ITRAN 


Source particle translation option ITRAN - 0, does 
not translate particle coordinates 
ITRAN = 1. translates particle coordinates by X , 

Y , (where X Q , Y q , Z q are entered separately) 
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TABLE X. - PRIMARY SOURCE GENERATOR INPUT 


Card 

group 

Format 

Card 

column 

Item 

Description 

a i 

7E9.4 

1 - 9 

10 - 18 
19 - 27 
28 - 36 
37 - 45 
46 - 54 
55 - 63 

E 

XO 

YO 

ZO 

ALPHA 

BETA 

GAMMA 

Energy in MeV of monoenergetic source 
particle when IE = 1; leave blank if 
IE * 1 

Source spatial coordinates translation 
vector. Leave blank if ITRAN = 0 

Monodirectional source particle direc- 
tion cosines when ED = 1; leave blank 
if ID * 1 

2 

(8E9. 4) 

1-72 

X(1,I) 

(1= 1,NX + 1) 

Values of FIRST spatial coordinate of 
the source 

3 

(8E9.4) 

1-72 

FX(1,I) 

(I = 1, NX + 1) 

Values of relative power distribution in 
FIRST spatial coordinate 

b 4 

(8E9.4) 

1-72 

FXS(1, 1) 

(I = 1, NX + 1) 

Values of relative biased power distri- 
bution in FIRST spatial coordinate 

5 

(8E9.4) 

1 - 72 

X(2,I) 

(I = 1, NY + 1) 

Values of SECOND spatial coordinate 
of source 

6 

(8E9.4) 

1-72 

FX(2, 1) 
(1= 1, NY + 1) 

Values of relative power distribution 
in SECOND spatial coordinate 

b 7 

(8E9.4) 

1 - 72 

FXS(2, 1) 
(1= 1, NY + 1) 

Values of relative biased power dis- 
tribution in SECOND spatial coordinate 

8 

(8E9.4) 

1 - 72 

X(3,I) 

. (1= 1,NZ + 1) 

Values of THIRD spatial coordinate of 
source 

9 

(8E9.4) 

1 - 72 

FX(3, 1) 
(1= 1,NZ + 1) 

Values of relative power distribution 
in THIRD coordinate 

b 10 

(8E9. 4) 

1-72 

FXS(3,I) 
(1= 1, NZ + 1) 

Values of relative biased power distri- 
bution in THIRD coordinate 

C 11 

(8E9. 4) 

1 - 72 

FAP(I) 

(1= 1, NANP) 

Fraction of source particles emitted in 
I** 1 polar angular interval 

C 12 

(8E9.4) 

1-72 

CP(I) 

(1= 1, NANP + 1) 

Cosine values of polar angles bounding 
angular intervals 

C 13 

(8E9. 4) 

1 - 72 

FAZ(I) 

(1= 1, NANZ) 

Fraction of source particles emitted in 
each azimuthal angular interval 

C 14 

(8E9. 4) 

1-72 

AZ(I) 

(1= 1, NANZ + 1) 

Azimuthal angle values (in radians) 
bounding azimuthal angular intervals 

d 15 

(8E9.4) 

1 - 72 

EG(I) 

(I - 1, NSGY + 1) 

Energy values bounding source energy 
spectrum 

d 16 

(8E9.4) 

1 - 72 

ECP(K) 

(K = 1, NSGY) 

Fraction of source particles emitted in 
each energy interval 

C 17 

(8E9.4) 

1 - 72 

EBP(K) 

(K = 1, NSGY) 

Biased fraction of source particles emit- 
ted in each energy interval 

i 


‘ l Omit card 1 unless IC = 0, or ID = 1, or IE = 1. 
d Omit unless IC 0. 

c Omit cards 11 to 14 unless ID - -2 or ID = 3. 
d Omit unless IE = 2. 
e Omit unless IE * -2. 


TABLE XI. - RESPONSE FUNCTION INPUT DATA 


Card 

group 

Format 

Card 

column 

Item 

Description 

1 

(12A6) 

1-72 

RSP(J, I) 
(I = 1, 12) 

Alpha-numeric identification for 
the response function 

2 

(8E9.4) 

1 - 72 

RESPFN(J,K) 
(K = 1, IEP) 

Response function data for the J th 

response at the K th output energy bin 


TABLE XII. - CAPTURE GAMMA-RAY SOURCE GENERATOR INPUT 


Card 

Format 

Card 

Item 

Description 

group 


column 



1 

715 

1-5 

NEC 

Total number of elements for which gamma 





rays are born from neutron capture 



6-10 

MNIC 

Maximum number of incident neutron energies 





in capture energy range for given element 



11 - 15 

MNGC 

Maximum number of different gamma -ray energies 





resulting from neutron captures for a given element 



16 - 20 

INEL 

Problem type: INEL = -1, inelastic gammas only; 





INEL = 0, capture gammas only; INEL = +1, 
both inelastic and capture gammas 



21 - 25 

IOP 

Option on number of capture gammas produced: 





IOP = 0, capture gammas are picked at random; 

IOP = 1, each neutron collision produces gamma ray 



26 - 30 

IETPE 

Tape number for inelastic gamma sources 

. 


31 - 35 

ICTPE 

Tape number for capture gamma sources 

a 2 

(319) 

1 - 9 

J 

Element number 



10 - 18 

Nil 

Number of incident neutron energies for element J 





over capture energy range 



19 - 27 

NGG 

Number of gamma energies for element J 


(3E9.4) 

28 - 36 

BLANK 

Leave blank 



37 - 45 

ELMAX(J) 

Maximum neutron energy (MeV) at which captures 





occur for element J 



46 - 54 

ELMIN(J) 

Minimum neutron energy (MeV) at which captures 





occur for element J 

a 3 

(8E9.4) 

1 - 72 

ENC(J, L) 

Incident neutron energy values (MeV) for which 




(L = 1, Nil) 

neutron capture data are given for element J 
(Loaded in ascending order) 

~ a r 

(8E9. 4) 

1 - 72 

EGAMC(J,K) 

Gamma energy values (MeV) resulting from 




(K - 1, NGG) 

neutron capture in element J (ascending order) 

a 5 

(8E9. 4) 

1 - 72 

GPC(J,L,K) 

Number of gammas, N. ; . , of k th energy released 
when neutron of incident energy L is captured by 





element J (Read in all gamma energies for the first 
incident neutron energy, followed by all gamma energies 
for the second incident neutron energy, and so on. ) 


a Cards 2 to 5 are read in for each element where capture gammas can occur. Reading 
order is cards 2, 3, 4, 5 for the first element, cards 2, 3, 4, 5 for the second element, 
and so on. 



Capture ga mma- ray source input. - The secondary gamma source generator re- 
quires a set of fixed point control parameters to determine the storage required for the 
secondary gamma problem. A table of secondary gamma production cross sections are 
also required. The format for preparing these data are shown in table XII. 


SUMMARY 

A modification and reprogramming of the COHORT general purpose Monte Carlo 
shielding program has resulted in a code that has the following features: 

1. The code can track either neutrons or gammas in any shield geometry that can be 
described by quadratic surfaces. 

2. Complete results from source generation to final output are obtainable in one 
single computer run for primary neutrons or gammas. For secondary gamma analysis, 
a source tape containing the secondary gamma source parameters is generated during 
primary neutron tracking for later re-use. 

3. Biasing features include the exponential transform during tracking and the selec- 
tion of primary source particles with respect to preferred locations, energies, and 
directions. All biasing is optional. 

4. Output is in the form of number fluxes sorted according to energy for volume 
detectors or number fluxes sorted according to angle and energy for point detectors. 
Standard deviations on all fluxes are computed. Quantities such as dose rate or heating 
rates are also obtainable from input response functions. 

5. The code is completely written in FORTRAN-IV. Variable dimensioning is em- 
ployed to make maximum use of computer memory. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, October 22, 1970, 

120-27. 
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APPENDIX - SAMPLE PROBLEM INPUT AND OUTPUT 


This section lists the input data and printed output from the code for a sample pro- 
blem. The sample problem is that of a point isotropic monoenergetic source of 3-MeV 
gamma rays in an aluminum sphere with a radius of 20 mean free paths. Point detec- 
tors are located at 1, 2, and 4 mean free paths. Region averaged fluxes are also com- 
puted. 

Program storage requirements are first printed out followed by a listing of the pro- 
blem control parameters. Next is printed a listing of the atomic weights, atom densities, 
the output energy, polar angle, and azimuthal angle bins, followed by a listing of the 
point detector coordinates. The sample output is shown under that heading. 

The library data are then printed out. Libraries 1 and 2, giving the boundary and 
region descriptions, respectively, are first printed. 

This is followed by a printout of library 7 which describes how the exponential trans- 
form is applied in each region. The exponential transform is not applied for the sample 
problem. 

Next is printed a listing of library 6 for each element present. This is followed by 
tabulations of internally computed probabilities and macroscopic cross sections. This 
printout may be deleted if desired. 

The control parameters for the primary particle source generator are next printed. 
The energy of the particles generated are listed. Since polar angle biasing is used in 
starting the photons for the sample problem, the listing gives the polar angle bins and the 
fractions of the particles assigned to each bin. A statement then follows that a source 
tape has been generated. The number of particles, total weight of the particles, number 
of records, and type of particle are listed. Where preferential selection of source 
particles is used, as for the sample problem, tables showing the distribution of the 
particle number and weight are then listed. In the sample problem, only polar angle 
biasing is used and a table of the polar angle bins and the number of particles in each bin 
is given. The azimuthal angle table shows that no source biasing has been used for the 
particles azimuthal angle. 

This completes the printout of the input data. The next two statements, BEGIN HTA 
and HTA COMPLETED FOR 31, indicate, respectively, that the tracking routines have 
been entered and have been successfully completed for the number of particles requested. 

The output data are then printed. First is a table listing how many particle histories 
have been terminated and for what reason. 

This is followed by a table giving the total energy deposition in each region in units of 
MeV per source particle. The energy deposition is broken down into energy deposition 
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due to particle absorption, scattering, or particle terminations due to minimum weight 
or energy. In addition, the total number of collisions in each region are also shown. In 
addition to the energy deposition, the standard deviation and the percent deviation (stan- 
dard deviation times 100 divided by the value) of the energy deposition are also given. 
Totals for each region and each process are also computed and listed. The same table 
is then repeated in units of MeV per cubic centimeter per source particle. This listing 
is valid only when the correct volume for a region has been entered, as in the case of 
regions 2, 4, and 6. Regions 2, 4, and 6 are thin spherical shells where region averaged 
fluxes correspond to those of detector points 1, 2, and 3, respectively. 

A listing of the output energy grid and any response functions entered is then given. 

Next are printed the region averaged fluxes for each region as a function of the 
energy bins selected. The flux, standard deviation, and percent deviation are given. 

The region averaged fluxes are then multiplied by the response functions present, 
and a tabulation of the region averaged response functions are then listed. 

Flux and response function output for the point detectors are then printed out. The 
detector number and its spatial coordinates are given. Fluxes at point detectors are 
sorted according to the energy bins previously read in and are also broken down into un- 
collided (UFLUX) and scattered flux (SFLUX) components. These are added to give a 
total flux (TFLUX). Standard deviations and percent deviations are also listed. 

A listing of the final random number used completes the problem output. The exe- 
cution time on the IBM 7094 for 1000 histories of this problem was 2. 9 minutes. 
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Sample Problem Listing 


POINT ISOTROPIC,! J l ASt 0 ) SJJRCE JF 
D3M31676 IOI 


3.3 M5y 3AHHAS IN ALJHINJH 


23 
d 1 

1 .0EO3 
J-3 
26. *>3 
.0633 
.? 

3. 

0.0 

0.0 


B 7 

12 3 

1.3 E - 3 6 


2 

1 

3.3 


25 

3.3 


3 l 3 
3.1 3.331 


103. 

360. 


RND NO 
SDN TL 1 
CDNTL 2 
C3NTL 3 
S3NTL 6 
AT* r 
UN 

2.339 SPRINT 1 
E 0 R IN T 
G»R I M T 
PPR INT 


23. 

3 * 

3. 

13.5 

20. 

3. 

3. 

21 . 

0 L T 1,2 

23. 

0. 

3. 62. 





-.1 ►. I 3 


1 T 







l. I «M l 


? 1. 

1 . 

L. 




13 3.3 t 

HD 1 


2 1. 

1. 

L. 




121.33 

in ? 


? 1. 

I . 

L. 




6 2 3.25 

9) \ 


2 1. 

1 . 

1. 




66?. > 5 

9 1 •, 


2 1. 

1. 

1. 




1 72 2 . > 5 

rtf) 5 


2 1 . 

1. 

L. 




U ) 6 . 2 5 

DO 6 


2 1. 

1. 

L. 




1.6 ♦ 

^8) 7 


2 8 







LIDS 2 

1 

1 1. 

1 2 






* E ^ 1 

2 

1 62.9 73 

2 3 

-1 1 





RES 2 

2 

l 1.3 

3 6 

-2 2 





RrG 3 

2 

1 131.967 

6 5 

-3 3 





R F S 6 

2 

1 1.0 

5 b 

-6 6 





R- 5 5 

2 

1 263.896 

5 7 

-S 5 





REG 6 

2 

1 1.3 

7 9 

-6 5 





R C J 7 

- 1 

f) 1.3 

-7 7 






r r j b 


6 26 

1 






LI0R 5 

15. 

.966 

.966 

.69 

13. 

1.326 

1.326 

. 653 

6 U 1 

0. 

1.3 78 

1.078 

. 773 

5. 

1.231 

1.231 

. 952 

5 AL 

5. 

1.166 

l. 166 

1.07? 

6. 

l. 3975 

1.3376 

l . ? 6 7 

6 AL 3 

3. 

1 .582? 

l. 5B22 

1 .6 y& 

2. 

1.9)3 

1.933 

1.933 

5 AL * 

1.5 

2 .2376 

2.2376 

2.2) 

L • 

2.75 

2.75 

?. 7> 

5 AL > 

.8 

3. 36 

3.05 

3. J5 

. 5 

3.6 3 

3.63 

3.63 

6 AL 6 

.5 

3.76 

3. 76 

3. 76 

. 6 

6. 1 3 

6 .12 

6.12 

6 AL 7 

. 3 

6.6? 

6.6 

6.6 

. 2 

5. 37 

5.29 

5.29 

5 AL 8 

. 15 

5.98 

5.77 

5. 7 7 

. 1 

7. 1 9 

6.61 

5. 61 

6 AL 9 

.03 

8.33 

6. 7? 

6. 72 

.36 

11.1 

7 . 1 

7.1 

6 AL 13 

. 35 

16.59 

7. 29 

7.2) 

.36 

22. 71 

7.51 

7.51 

6 AL 11 

.33 

66. 76 

7. 76 

7. 76 

.02 

169.3) 

3.33 

8.33 

5 AL 12 

.315 

353. IH 

8 . 1 8 

8 . 1 8 

.31 

1 18 3. 

3.3? 

9. 32 

6 AL 13 


7 H 







L18R 7 


1 3. 

L 

3. 


1 3. 

1 

3. 

EXPIRE 1 


1 3. 

1 

3. 


1 3. 

1 

3. 

EX n TRF 2 

3 

-2 1 

J 3 

3 3 

-6 

1 j 



sdj <; c 

3. 








t 

. 1 

. 35 

. 25 

.15 


l .35 



PDL FRCT 

1 . 

. ->998656 

.996196 . 

965933 

. 46t»0 3 

. 7371 

-1 . 


CDS T H ( T 

1 . 








A 7 H FRCT 

3. 

6.28 







A 7 M AN(,L 


HA^HA OISE 

R A TF ( R/ HR PER PH )TjN/_ 

hp-se: j 




.36 - 

J 1.3526-6 

l . 6 l 2 9- b ; 

'. 2 72 7-6 

P.8 795- 

S 3.5 33-5 

6.38 16-5 

6 . 3 3 5 - 

6R/HR 


3 3 3 3 ) f L= 3)332 


Sample Problem Output 


p MNT I S ;T-O d I C t HI a bEJ) SJJR15 3^ ).> M Zv In AlJMlNJH 


9.5 * STHAOt RtOJIKCO 'JR PRDJ.EH =k”JM3n 
12*3 * STDRAGh RtUJIRED DJT*j! 


IMIIAI 
I DP 
t-PS/ 
MT 
M80 

r-iLii-i 
OX SEC 
JTLF 
I E P 
I PP 

I S T A J 
Kb I A$ 
•SK, 
HIST 
t M I N 

ElO 

HR 


R N3 l v ND. =033 1)16 7<*10 1 
PR id . . y I) « j. 

ND. 3- p AP T I C I • i 1 DOP 
NO. 3- H A T c R I ft i-S 
ND. J- rtO JN -ARIES 


PR J8LEH S 3 N T R )L S 


ND. 3 I N° J T L 10RAR It S 
CR jS> -^FCTI i'* •»-»! NT OPTION 
Rf .. : L J X F > j/1/2/3, NJ/V-S/1HAI. /, 

\J. J : IjTPJI tNERSY ft!\S 
NJ. 1- )JTj jt A? I H J T hAl. AmD! E 31 N S 
SJtRE- TAPE OPT I ; H 3/K, INTERNA,./? A-* 
FxrD. TRANS, hi 1ST FL T . 0 / l, AL l / 1 V 
CjILtilDN Tape OPllDN 3/1 
N*J. 3 - ilSTD-UfcS 
MM. / JT DEE ENt R,,Y< HTy ) 

EN 1 - R jf OfcPjSITIDN SJT'JFEIhEV) 

R J > S l AN ODJLETTE P4 RA«fTE*» 

DTI Mr CJT 3= : T I He I Hi N. I 


l/l , N 3/ Y: 


N3/ 1- S 


2 D 

nr ; 

pari i : . ; 

-> t 

NrL 

N t. )= . 

L 

♦ RED 

N - 

7 

N) )MX 

«U. NJ. 

6 

N ) Y 

N j , J' 

:S 1 

N ) 

N t. ) : ■ 

MM3- 1 

LIND 

L ■ A < A ) f 

9 

l S J 

N). J = 

1 

NS D R 

nj. d ; : 

: < D 

MAXCjL 

H A X . ND. 

T DN„Y D 

D 

I S 3 3N 

NJ. 3 1 ■ 

1.3303503 

f M A X 

H A X . E N i 

2. J003E-31 

wHI N 

PARTI )L: 

3. 

t T M 

T Hr RHAL 

1.33335-31 

7 »SL 

* p S I L 3 N 

3. 33D3F- J1 

tl: jt 

ND. H.F, 


I Yi»- jt l, N/ 

*ENTS 
.1 'NS 

JE 3 JJN3AR I Ei / Rc D! 'N 
:NE R .l : i IN M. 6 
J 3 I \ T J - TMI JRS 
TAPE N 3. 

J JTPJT ’ )L A R AnS.E 3 INS 


WT. SJT3FP 


^3 ^EH-DO 
i -36 


3E3H. CA.IS 


l- 3 3 
♦ 01 



ATWT(NEL) 

26. 993 

UNI MAT.NEL) 

6.0300E-32 

ESIINEL) 

0 . 

OUTPUT ENERGY BOJNDS 

2.030E-31 4.03DE-QI 6.000E-31 1.330EO3 1.500603 2.333EOD 2.5D3E+33 2.999EO0 3. 333EOO 

output polar angle bojnds 

0. 183.000 

OUTPUT AZUJTHAL ANGLE B3UNDS 
0. 360.000 


point oete: r or locations 


L 

fp: i l ) 

XDK. ) 

YD 11 L ) 

ZOllL) 




l 

2.333E+01 3. 

0. 

1 • 053E+01 




2 

2.333E+01 0. 

0. 

2 . 103E + Q 1 




3 

2.333E+31 0. 

0. 

4.200E *-01 




LIBRARY 

TY 3 E 

1 boundary description 





NO. 

IBM J) 

AF ( J ) 

BFt J) 

:f ( j ) 

XF C J 1 

Y F < J 1 

ZF(J) 3F ( J ) 

1 

2 

1 . 0030E+00 

1 • 0300E+03 l . 0003E + 0 3 

“3. 

-3. -3 

1 . 3303 E +3 2 

2 

7 

1.003 3E+00 

1. 000QE + 03 1 . 0000E+03 

-0. 

-3. -3 

I. 21 03 E +02 

3 

7 

1.0330E+00 

1.0300E+00 1 • OOQOE +03 

-3. 

-3. -3 

4.2325E+02 

A 

2 

1.0030E+00 

1 . OOOOE <-00 1. 3033E + 30 

-3. 

-3. -3 

4.62256*02 

5 

7 

1.0033E+00 

1. 3000E+03 1.3033EO3 

-3. 

-3. -3 

1.72226+03 

6 

? 

1 . 0D33E+00 

l.OOOOE+OO l . 00Q0E+03 

-3. 

-3. -3 

1 . 83 62 E +0 3 

7 

2 

1.0030E+00 

l.OOOOE+OO 1 . 0 OOOE + 33 

-0. 

-3. -3 

1.63D3E+35 

LI BRARY 

TY 3 = 

2 REGION 3ESCR IPT I ON 





NB 

mt 

VOL. I I bn 

MP1 





1 1 

L 

1 . OODOE+OO 1 

2 





2 2 

L 

6.2973E+01 2 

3 -1 

1 




3 2 

L 

1 . 3000E *00 3 

4 -2 

2 




A 2 

l 

1.3195E+02 A 

5 -3 

3 




5 2 

L 

1 . OOOOE *00 5 

6 -4 

4 




6 2 

L 

2 • 6389E+02 6 

7 -5 

5 




7 2 

l 

1.3003E+00 7 

8 -5 

6 




8 -1 

3 

l.OOOOE+OO -7 

7 





LIBRARY 

T Y 3 E 

7 EXPONENTIAL TRANSFORM 





REGIONS 

1 

2 3 

4 5 

6 7 

B 



NEXT! J) 

1 

l 1 

1 1 

i 1 

1 



EXT(J) 

3. 

0. o. 

0. 3. 

0. 3. 

3. 



LIBRARY 

TY 3 E 

6 CROSS-SECTIONS 









I 

ENERGY! I 1 

r n i j » i ) 

SM( J * I ) 

ESMl J, I) 

ELEMENT 1 

N0.= 1 

ATWT(J)= 26.983 









1 

3. L50003Ef02 

D.964330EO3 

3.9S4333E03 

0 . 49303DE+03 




2 

0. 103003E+32 

3. 1 32 433 = + 3 1 

3. 132433E01 

3.663 ODE+33 




3 

3. B33303EOI 

3.137333E+31 

3. 1 379 3 3 E + 3 l 

3.773033 t+33 




4 

0.63000QE+01 

3.128103EO1 

3 . 12B1 33E + 31 

3 . 95 2 U33 6+03 




5 

3. 533300E + 31 

D.116433E+31 

3. 116433E+31 

3.1D7733E+31 




6 

0. 433003 = 01 

3.138763EOL 

3. 13B7S3E01 

3.124 703 E + 0 1 




7 

3. 300000E+ 31 

3. 1 5822 3E + 3 1 

3. 158223EOL 

3.149633E+31 




8 

0. Z00003E+D1 

0 . 1 93 3 DDE + 3 l 

3. 1 9333 3 E *-3 1 

3. 193 33 3 E+3 1 




9 

0. 153003E+01 

3. 223763=+3l 

3. 2237S3E01 

3.223033E+31 




13 

3. lOOOOOEf 01 

3.275330E+31 

3. 275D33EO! 

3.275030E+31 




11 

3. 830000600 

0. 30630DE + 3 l 

3. 306D33EO1 

3*336 003E+31 




12 

0.600003E+0D 

0.3483035+31 

3.348333E+31 

D.34B033E+31 




13 

3. 530003 = 03 

0. 376333 = 01 

3. 376333E+31 

3.37603DE*31 




14 

0. 400000E + 03 

3.413330EOL 

3. 41 2 3 33E + 3 1 

3.4120336+31 




15 

3. 3000006+33 

0. 462300Ef 31 

0.4S0033E+31 

D.463033E+31 




16 

0.2D0003E+33 

0.5373335+3 L 

3. 529333 = 01 

3.529033E+31 




17 

3. 153003E+D3 

3. 598333=01 

3. 577333E+31 

3. 577033 £+31 




18 

3. 1 0 OOOQE + 33 

3.719333E+3L 

3.641033E+31 

3.6410D3E+3 1 




19 

3. 800000E-31 

3. B330335 + 3 l 

3.672033E+31 

3.67203DE+31 




20 

0-630000E-31 

3. 1 l l 3335+32 

3.713333E+31 

3.713033E+D1 




21 

3. 5330D3E- 3 1 

3. 1 45 93 3E + 3 2 

3. 729333E+31 

3.729033E+31 




22 

3. 40D300E-01 

3.2271 33E+ 32 

3.751333E+31 

3.751033F+31 




23 

3. 300000E-31 

0. 457603E+32 

3.776333E+31 

3.776033E*3l 




24 

0. 2D3333E-31 

3.149333E+33 

3.33333 3E01 

3.933 03 3E+31 




25 

3. 153303E-01 

3. 353183 = 03 

3. 318333E+31 

3 .313033 E+31 




26 

3. 130000E-31 

3.118033EO4 

3.832333E+31 

3 . B32 333 E* 3 1 
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COLLISION PROBABILITY TABLES 


ENERGY(NGY), NGY = 26 


1.5000E+DI 1 . 0 OOOE + 0 l 

8 . OOOOE-O 1 S.OOODE-Dl 
5. 0000E-D2 *• 0000E-02 

3.00006+00 
5 • OOOOE-O 1 
3.0000E-02 

6. OOOOE+DO 
*. OOOOE-Ol 
2. OOOOE-O 2 

5 • 330 0= H-00 
3.0000E-01 
1.5000E-02 

*. ODDDE + D D 
2. 0DDD5-DI 
1.0D00E-D2 

3.0D00E+DD 
l. 5300E-DI 

2-DDOOE+OO 
1. DDDDE-D1 

1 • 5DOOE+OD 
8.000DE-O2 

COMPTON/ ELASTIC SCATTER 

0- 508299 0.6*7*61 

1.000000 1.000000 

1- 000000 1.000000 

PROBABILITIES 

0.721707 

1.000000 

1.000000 

FOR ELEMENT 

0. 7*3159 

1. DOOOOO 
l.DOOOOD 

1 

0.925258 
i. 000000 
l.OOOOOD 

0. 39B57* 

1. DDDDDD 
1. OOODDD 

0.9*5519 
t. DDDDDD 

3.93**80 
I. 0DDD30 

0. 996503 

1. DDOOOj 

TOTAL MACROSCOPIC CROSS 
0.05813 0.06175 

0.18*52 3.209B* 

0.87973 1.369*1 

SECTIONS FOR 
0.05500 
0. 22673 
2. 81963 

MATL l 

0.0772* 
0. 2*93* 
8. 98651 

D. 07019 
0.27359 
21-29675 

D.DB357 
D. 32331 
71.15*00 

0.095*1 
D. 35059 

D. 11556 
D. *3356 

D. 13*93 
D. 50230 


PRIMARY SOURCE GENERATOR C3NTR3LS 


IC 

SOUR; r 

GEOM. D/l/2'3. 

PT./RECT./CYL./SPH. 

0 

IE 

SOURCE 

ENERGY 

DISTR IB. 

1/2 ,M3ND/SPECT 

1 

NY 

N3. D = 

SPATIAL 

INTERVALS IN 2ND DIR. 

0 

NSSY 

NO. D = 

ENERGY 

BINS IN 

SOURCE SPECTRUM 

0 

NANZ 

NO. 0 = 

AZIMUTHAL ANG. 

BINS IN S3URCE DIST. 

l 


ID SOJRCE ANG. DISTRI3. 1/2/3 » MONO/ I S3/ AN ISO 

NX N3. 3 : SPATIAl. INTERVALS IN 1ST Dll. 

NZ NO. 0= SPATIAL INTERVALS IN 310 DM. 

NANP NO. 0 C POLAR ANGLE BINS IN SOJRCE DIST. 
IRAN SOJRCE COORD. TRANSLATION 0/1 .NO/YES 


ENERGY OF M 0 N OEN ER GET I C SOJRCE 3. OOOQE+DOME V 


SOURCE FRACTION IN tACH POLAR INTERVAL 

l.OuDDE-Ol 3.5000E-01 2.5000E-01 


1.5000E-01 


1.0D00E-D1 5. 0030E-32 


COSINES OF AN3JLAR INTERVALS 

1.00005+00 9.9986E‘01 9.962DE-Q1 


9.6593E-01 


8.5603E-01 7.07L05-01 -1.DDD0E+DD 


SOURCE FRACTION IN EACH AZIMUTHAL INTERVAl 
1.00005+00 


AZIMUTHAL ANGLES 

0. 6 • 28Q0E+00 


SOURCE PARTICLE PARAMETERS HAVE BEEN GENERATED AND WRI TTEN ON SCRATCH T A 3 E NO. 2. 
THIS TAPE CONTAINS * 

I .00000 OE + 03 = NUMBER OF SOJRCE PARTICLES 
9.99*88DE+D2 = THE TOTA. WEIGHT OF THE SOJRCE PARTICLES 

11 = NUMBER OF LOGICAL RECORDS {MAX. OF 100PARTICLES PER RECORD) 

100 * NUMBER OF PARTICLES ON LAST RECORD 

l = TYPE OF SOURCE PARTICLE/ 0=N£JTRDNS/ 1=GAMMAS/ 


PO.AR ANGLE AND WEIGHT DISTRIBUTION FOR A BIASED ISOTROPIC SOURCE 


INTERVAL NO. 

C3SIN5 

THETA 

NO. OF PARTICLES 

ANGJLAR WEIGHT 

1 

1.0000 - 

0. 9999 

10D 

7.22 *3 2 3E-D2 

2 

0.9999 - 

0. 9962 

350 

1.3297*9E+00 

3 

0.9962 - 

0.9659 

250 

1.5131 *3 C +01 

* 

0.9659 - 

0. 0660 

150 

*. 9951 *3E+D1 

5 

0.8660 - 

0. 7071 

too 

7.9*5*96E+01 

6 

0.7071 - 

-1.0000 

50 

3. 535500E+02 



TDTALS 

1000 

9. 999999E +02 


AZIMJTHAL ANGLE AND WEIGHT DISTRI3UTIDN “DR A BIASED ISOTROPIC SOJRCE 


NTERVAL ND. 

ANGLE LI 

[ MI TS( RADI ANS ) 

ND. 3= PART1C.ES 

4NGJ.AR WEIGHT 

1 

0. 

6.2800 

1300 

9.99*3905+02 



TDTALS 

1000 

9. 99*390E+02 


BEGIN HT A 


.ODOOE+OO 
.00035 - 02 


l.OOD jOO 

1.000000 


0. 16582 
0.66933 


-2 

0 

0 

-6 

0 
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particle history terminations 

NUMBER OF 

PARTICLES REASONS FDR TERMINATIONS 


3 . 

L.0D0030EO3 

0 . 

0 . 

0 . 

0 . 

Z. 00003 JE+33 

4. 3310035+33 
4.3310DDE+D3 


NUMBER DF COLLISIONS = 12. THE MAXIMJM ALL3J53 

THE SCATTERED ENERGY IS LESS THAN 2.33035-31 M=V. THE MINIMJM AlLONED 
THE WEIGHT IS LESS THAN 1.0030E-Q5, THE MINIMJM ALLOWED 
ESCAPES FROM SOJRCE REGION JNC0LLI3E3 

ESCAPES FROM SOJRCE REGION AFTER DNE DR MORE COLLISIONS 
ESCAPES FROM SHIELD GEOMETRY 

THE NUM3ER DF HISTORIES ACTUALLY PROCESSED IN THIS PR3JLEM WITHIN THE CJT-OFF TIME LIMIT 
THE NUM3ER OF EXPONENTIAL CJT-DFFS IN TRAC< LENGTH ANA.YSIS 
TOTAL NJMBER DF COLLISIONS 


TOTAL FnERGY DEPOSITION BY REGION IN JNITS DF (MEV/SOJRCE ^ARTICLE) 


REGION NJMBER 


PROCESS 


1 




ENERGY 

STD OEV 

PC DEV 

cNERGY 

NO. COLS. 


1497 



ABSORPTIO N 

1.3256-33 

2.29IE-04 

2.234EOI 

B.658E-35 

SCATTERING 

1.4535+30 

1.944E-01 

1.32/5+31 

2. 737E-32 

MIN WT+ENGY 

4 . 9 3 3E-0 2 

I . OD 8£- 02 

2.043E+01 

6. 3795-33 

totals 

1 . 5 1 5E + 00 

2.002E-01 

1.321E+01 

J.354E-D2 




REGION 

NJMBER 

PROCESS 


4 




ENERGY 

STD DEM 

PC DEV 

ENERGY 

NO. COLS. 


87 



ABSORPTION 

1.544E-04 

9. 3776-35 

6.071E+01 

1.235E-33 

SCATTERING 

3.379E-02 

1.6376-32 

5.3l3E«-3l 

4.451E-31 

MIN WT+ENGY 

5.933E-03 

3.7586-33 

5. 3 3 4E + 0 1 

5.237E-32 

TOTALS 

3. 588 E- 02 

I .720E-32 

4. 655E+0L 

4.9B4E-D1 




REGION 

number 

PROCESS 


7 




ENERGY 

STD DEV 

PC DEV 

ENERGY 

NO. COLS. 
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ABSO^PTIIJ J 

7.332E-05 

4.2596-05 

5.9396 *3 I 

3. 

SCATTFRI N 

4. 165E-32 

2. 3106-32 

5.5456+01 

J. 

MIN * T + E N , Y 

L.B37t-03 

4. 2746 -34 

2. 3255+31 

3. 

TOTALS 

4. 356E-02 

2-3396-32 

5. 3316+31 

J. 



TOTAL ENERGY DEPOSITION BY 

REGION IN 

JNITS OF ( Mr V/ 




Rc G I J N 

N JM3ER 

PRUCf SS 


1 




energy 

STD DEV 

3 : dev 

ENERGY 

NO. COLS. 


1497 



AbS OR p T 10 • 

1.325E-33 

2-29 IE-34 

2.2346+31 

1.3765-35 

SCATTER IN , 

l. 4656 +30 

1. 9445-31 

1. 3276+31 

4.299>=-04 

M[ N WT+ENG 

Y 4.933E-32 

1 .038F-32 

2.343E+31 

1.31 -34 

TOT Al S 

1. 5I5r+33 

2.032E-31 

1.3215+31 

5. 3265-34 




REGION 

NJMBER 

PRCCt SS 


H 




ENERGY 

STD DtV 

PC DEV 

ENERGY 

NO. C°LS. 


87 



AriS.l *PT I D i 

1 . 1 7 j 6 -36 

7.1365-37 

5.071E+31 

1. 2355-03 

SCAi TtRlN • 

2. 3 3 3 c - 34 

L .2415-34 

5. 31 H + 'Jl 

4.4515-31 


2 



3 


STD DEV 

PC DEV 

ENERGY 

STO 3 5 V 

PC OE V 

161 



1254 


3.4595-35 

3.9B9E+DI 

1.331F-33 

2.5316-34 

1.977. +31 

9.1235-33 

3.3715+31 

3.113 c-3 1 

1. 732E-01 

2. 135*; + 31 

4.3515-33 

5.3515+Dl 

5. 6776-32 

1.3276-02 

2. 337: +J1 

1.339r-02 

3.932E+31 

8.6966-31 

1.8496-01 

2. 1276+31 

5 



6 


STD DEV 

PC DEV 

ENERGY 

STO DEV 

PC DEV 

1 33D 



31 


3. 3935-34 

2.5355+31 

5.3446-36 

2.35/6-36 

4. 377i +31 

1 • 42 Db- 31 

3. 191E+31 

i. 2566-03 

5.1 99E-D4 

4. 137- +31 

1. 441=02 

2 . 75BE+31 

2.8256-34 

1.312* - j4 

4.6425+31 

1. 535E-31 

3.3B3E+31 

1.5446-33 

5.9346-34 

3. / 7 B : + 3 1 

3 



T3TAL5 


STD 35V 

PC DEV 

ENERGY 

STD DtrV 

PC DEV 

D 

3. 

3. 

3.9136-33 

5.8976-04 

1. 5085+31 

3. 

3. 

2.922 +33 

3.167t-Ol 

1. 1225 + D1 

3 . 

3. 

1 • 7 2 6 3 1 

?.23?r-D2 

1.279E+31 

3. 

3. 

2.99BE+30 

3.175 t-:> 1 

1.3595+31 


1++3/SJJRC5 

PARTJC-5) 




2 



3 


STD 35V 

PC DEV 

ENERGY 

STD 05V 

PC DE V 

161 



• 264 


5«49l:-j7 

3.9395+01 

1.3316-33 

2.531E->4 

1. 977:+31 

1.453:04 

3.3715+31 

3.115 -31 

1. 732c-3 1 

2. 1355+j 1 

5.433-35 

5. 351E+31 

5.577 r -02 

1.327F- 2 

2. 3 576 + 01 

2.073- -34 

3. 9325+31 

3.596^-31 

l .94 ), - ] 

2. 1275+ 31 

5 



6 


STO 3 : V 

PC DEV 

ENERGY 

STD DEV 

PC DEV 

1333 



31 


3. 393: -34 

2.5355+31 

1.9UE-38 

7. 793E-39 

4.077 +31 

1.423 - 0 1 

3.191501 

4,7615-36 

1. 97jF-3s 

4. 1 37 + J1 
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I I I 


■ I I ■■ I ■ ■ ■■ I 


min wt*engy 

4.496E-05 

2. 8485-05 

6 . 3 3 4E *0 l 

5.237E-02 

1.4416-02 

2.758E*0l 

1.071E-36 

4. 971 t -07 

4. 6426*01 

TOTALS 

2.7956-04 

1 .304E-04 

4. 665E* 01 

4.9B4E-01 

l. 5366-01 

3. OB36*01 

5. 951 F-06 

2.211E-06 

3. 779i*0l 

PROCESS 

ENERGY 

7 

STD OEY 

REGION 
PC DEV 

N JMBEiR 

ENERGY 

9 

STO 3 = Y 

PC 05/ 

ENERGY 

TOTALS 
STO DEV 

PC 3 = V 

NO. COLS. 
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3 





ABSORPTION 

7.332E-05 

4.259E-05 

5. B 09E ♦ 01 

3. 

3. 

D • 

3.666 E- 33 

5- 897E-04 

1. 5086*01 

SCATTERING 

4. 165E-02 

2.310E-02 

5.5466*01 

3. 

0. 

3. 

2. 754EOO 

3. 167E-31 

1. 122E+31 

MIN WT+ENGY 

I.B37E-03 

4.274E-04 

2. 32SE01 

3. 

3. 

0. 

1.602E-3 1 

2. 237L-02 

1.279E+01 

totals 

4.356E-02 

2.309E-02 

5.301Ef 31 

3. 

3. 

3. 

2.927 EOO 

3. 1 75E-01 

1.0596 + 01 


response function inpjt 

OUTPUT ENFR37 BOUNDS 


2.0000E-31 4. 0000E-01 

5.0000E-31 1.0030E*30 1.5000E*00 

2.3D33EOO 

2. 5003E+33 

2.9990E*00 3.0030E+00 

GA«HA DOSE RATE 
6 . 060E-0 7 1.353E-06 

(R/HR PER PHOTON/:N2-SEC ) 
1.613E-06 2.273E-06 2.8996-05 

3. 539E-3S 

4.082E-05 

4. 336E-06 



FLUXES VERSJS REGION ANO ENERGY 
LCHER ENERGY = 2.000E-01 ME V 







UPPER ENEGY 
BOS ( MEV) 

FLUX 

STD DEV 

REGION 
PC 3EY 

NJM9 

ER 

flux 

2 

STD DEY 

PC DEY 

FL JX 

3 

STO DEV 

RC OEY 

A.ODOE-jI 

1.5516*33 

3.465E-31 

2.23AE*01 


2.613E-33 

1.3B2E-33 

4.140E*3L 

1.9796+30 

4, 459E-31 

2.2541*01 

5. OOOE-Ol 

1 • 3 9fa E *00 

2. 76 AS -31 

2. 52?E*31 


2.343E-33 

1. 033E-33 

4.28BE*31 

l. 191E+00 

2.537E-01 

2.2146*01 

1.300E* JO 

4.544E-31 

l .936E-31 

4. 1 9 5E * 3i 


9. 1115-34 

4 • 293 c ~0 4 

4.709EO1 

7.D42E-31 

2 * 6 31 F -01 

3.7356*31 

1.5006*30 

4. 31BE-31 

1.771E-01 

4.4086*01 


1. 4395-03 

6. 921E-34 

4. SIDE *31 

7.087E-31 

3.554E-01 

5.3 156+01 

2 . 0 00 E *0 3 

1.872E-31 

9.779E-32 

5.225E+D1 


3.451E-34 

1. 7196-34 

4. 9836+01 

2.696E-01 

1.072E-01 

3. 9756*01 

2.5006*03 

5.512E-31 

2.339E-D1 

4. 1 39E* 31 


9.121E-34 

4. ID9E-04 

4.5346+01 

3.948 E- 01 

1.4556-01 

3. 6856*01 

2. 9996+33 

2.5166-01 

1 .A89S-01 

5.92DE*31 


8. 336E-04 

5. 944E-34 

7.3966+01 

3. 0946-31 

2. 191E-0 L 

7.384E+31 

3-300E+30 

5.4416*00 

7.265E-31 

1.128E*01 


5. 833E-03 

6. 576 E -04 

1.1296+01 

2.196E*33 

2 • 4 6 5 t -01 

1 . 1286*01 

TOTALS 

l . 3 93E *31 

1 . 295E + 30 

1. 184E*01 


1.5196-02 

2. 543E -33 

l. 73 8 E * 3 1 

7. 742fc*30 

1. 286E+0 

1.661 E* 01 

UPPER ENEGY 
BDS l MEV J 

FL JX 

A 

STD DEY 

REGION 
PC OEY 

NUMBER 

FLUX 

5 

STO OEY 

PC DEY 

FL JX 

6 

STO DEV 

=>C 3 c Y 

4.030E-31 

1.331E-33 

3. 48 IE-34 

3. 473E*01 


1.8336*33 

5. 239E-31 

2.7556*31 

1. 264E-34 

9. 953E-35 

7 . 869 E* 0 1 

6.330E- jl 

6. Z 84E-04 

1.B13E-04 

2. B35E+01 


1. 358£*33 

Z.595E-31 

2.547E+01 

2. 698 E- 35 

9.012 E-OS 

3. 3406*01 

1.0006*03 

2.984E-04 

1 . 333E-OA 

4.366E*31 


3.9366-31 

1. 327E-01 

3.372E+31 

l. 9416- 35 

5. 5 72E-06 

3. 3866*01 

1.5006*33 

3. 33BE-04 

1 -386E-3A 

4.153E*3l 


2.743E-31 

9. 5056-32 

3.532E+01 

1.313E-35 

3.452E-0‘. 

2. 534:*31 

2.333E+03 

1.9126-04 

9 . 395E-05 

4. 914E*01 


3.993E-31 

1. 5946-01 

3.9916+01 

3.D24C-D5 

1.249E-05 

4. 1306*01 

2.5006*33 

2.175E-04 

8.838E-35 

4. 064E* 01 


4. 885E-D1 

2. 333E-31 

4.7B5E+01 

3. 409 t- 35 

1.7776-33 

5.2156*01 

2.9996*33 

1.471E-04 

9. 9466-05 

6.7636*01 


1 . 743E-31 

1 • 080E- 01 

S. 1986+31 

1.025E-35 

6 • 1 60 t -06 

6.0166*01 

3-0006*33 

l. J22E-33 

1.153E-34 

1. 128E*01 


1 . 1 47E*03 

1. 294 E- 01 

1.1286+01 

5. B93E-35 

7.771 E -05 

1.1 286*01 

TOTALS 

3 . B39E-33 

6.165E-34 

1. 636E*31 


5.813E*30 

1.2336*30 

2.1286+01 

3.294 E- 34 

1 . 102E -34 

3. 3 44 E* 0 1 

UPPER ENEGY 
RDS(MEV) 

FL JX 

x 

STD DEY 

REGION 
RC OEY 

NUMBER 

FL JX 

8 

STO OEY 

PC DEY 




A.390E-J1 

l. 374E-01 

1 -2846-01 

9. 347E*0I 


0. 

3. 

0. 




6 • OOOE-? 1 

1.701E-03 

l . 701E-03 

1. 333E*02 


3. 

0. 

3. 




1 . OOOE+33 

0. 

3. 

3. 


3. 

3. 

3. 




1.5006*33 

3. 

0. 

3. 


3. 

0. 

3. 




2.000E*03 

3. 

3. 

3. 


3. 

3. 

3. 




2.500E+33 

0. 

0. 

3. 


3. 

3. 

3. 




2.999E* ) 3 

0. 

3. 

3. 


3. 

3* 

3. 




3.O0DE* O 

3. 

3. 

3. 


3. 

3. 

3. 




totals 

1.391E-31 

1 . 284E-01 

9.2286*01 


3. 

3. 

3. 
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GAMMA DOSE RA TE { R/ HR PER PHOT 3N/ C M 2-S5 C ) 
LOWER ENERGY = 2.000E-01 

RESIDE NUMBER 


UPPER ENEGY 
BQS(MFV) 

FLUX 

1 

STO DEV 

p: 3EV 

FLUX 

2 

STD 3EV 

PC DEV 

FLJX 

3 

STO DEV 

PC 3 _ V 

4 . 030E-0 i 

9 . 397E-D7 

2. 130E-37 

2.234E+D1 

1.5B4E-39 

6.5555-13 

4. 1435*01 

1.199E-36 

2. 732)—37 

2.254:01 

6 • DOQE-Ol 

1.154E-06 

2.909E-37 

2.522E+31 

2.463E-39 

X.355E-39 

4.28B501 

1.254E-36 

2.776E-07 

2. 214 E*- 31 

i.oooe+oo 

7.329E-07 

3-075E-07 

4.195E+01 

1.473E-39 

6.923E-X3 

4. 7O9E+01 

X.135E-36 

4.243E-3 i 

3.7355*31 

1.500E+D3 

9. 132E-37 

4 . 025E-07 

4.43BE+Q1 

3.270E-09 

1.5735-39 

4.813E*-31 

1.51XE-D6 

3.3736-07 

5.C15;+31 

2.000E+00 

5 . 425E-3 7 

2.835E-37 

5.225E+01 

1.333E-39 

4.9315-13 

4. 9B05+01 

7. 815 E- 37 

3.13SE-07 

3. 9755*31 

2.500E+00 

1. 934E-06 

B.131E-37 

4.139E+31 

3. 233E-39 

X.442E-39 

4. 504E+31 

X.385E-36 

5.135E-07 

3. 6855*01 

2.999^+3 3 

1.327E-36 

6.O78E-07 

5.923E+01 

3.2B3E-39 

2.425E-39 

7.395E+31 

1.263 E- 36 

3.945T-07 

7.3845+31 

3.300E+33 

2.325E-D5 

3.187E-06 

l. 1 2BE+01 

2.557E-0B 

2.8B4E-39 

1 • 1 2BE+31 

9.537E-36 

1.381 E-36 

1.1285+31 

TOTALS 

3.549E-05 

4. 20 2E-D6 

1. 1 B 9- E ** D 1 

4. IB4E-33 

7.2735-39 

1. 738EOL 

i .B22E-35 

3.327E-06 

1.6615*31 

UPPER ENEGY 
BDSIMEV) 

FLUX 

4 

STD DEV 

REGION 

p: 3E V 

NUMBER 

FLJX 

5 

STD 3EV 

PC DEV 

FLJX 

6 

STD DEV 

PC 35V 

4. OOOE-31 

6.364E-10 

2.109E-10 

3.47SE+31 

1. 141E-36 

3.157E-37 

2. 7665+31 

7.662E-11 

6.3336-11 

7. 869E+31 

6. OOOE-31 

6.614E-10 

1.908E-10 

2-8B5E+01 

1.114E-36 

2. B36E-37 

2.547E+31 

2.840E-11 

9.486E-12 

3. 343E+D1 

1. 0006*33 

4.8I2E-10 

2.X01E-10 

4. 3 66E + 01 

6.348E-07 

2.143E-37 

3.372E+31 

3. 131 E- 11 

1.363E-11 

3.3865+31 

1.530E+33 

7.5B7E-10 

3.149E-10 

4.153E+01 

6.233E-07 

2. 133E-37 

3.5325+01 

2.978 E- 11 

7.345E-12 

2.6345+31 

2 . 000E+03 

5.542E-10 

2.723E-X0 

4 • 9 1 4E +0 1 

1.157E-36 

4.519E-37 

3.991E+31 

8. 766 E- 11 

3. 523E-11 

4. 1335+J 1 

2.500E+03 

7.631E-10 

3.X01E-X0 

4.064E+01 

1. 7146-36 

8. 232E-37 

4. 7 85E*-31 

1.196E-13 

6.2376-1 1 

5.2155+31 

2.999E+33 

6.306E-10 

4.06DE-10 

6. 76 3E * 0 1 

7.115 E ~07 

4. 41DE-37 

5.193E*Ol 

4. 135 E- 11 

2.5176-11 

6. 31 5- + 3 l 

3.030E+33 

4.4B2E-09 

5.055E-ID 

1.12BE+01 

5.031E-36 

5.6745-37 

1.1235+01 

3.3226-10 

3.4336-11 

1. 123: + 01 

TOTALS 

B -937E-09 

X * 430E-39 

l. 505E +31 

X. 213E-35 

2. 5835-36 

2. 128E+D1 

7.174E-1 3 

2.399E-1 

3. 3445+31 

UPPER ENEGY 
BDSIMEV) 

FLUX 

7 

STD DEV 

REGION 
PC 3EV 

NJMBER 

FLJX 

3 

ST3 0 5V 

PC DEV 




4 .OOOE-31 

8.328E-08 

7.784E-08 

9.347E+31 

0. 

0. 

3. 




6. OOOE-31 

1. 790E-09 

1.790E-09 

1.333E+02 

3. 

3. 

3. 




1.000E+03 

0. 

0. 

3. 

3. 

3. 

3. 




1.500E+03 

0. 

0. 

3. 

3. 

3. 

3. 




2.000E+00 

0. 

0. 

3. 

0. 

0. 

3. 




2.500E+33 

0. 

0. 

3. 

3. 

3. 

3. 




2 . 999E+03 

0. 

3. 

0. 

3. 

3. 

3. 




3 • OOOE+ 3 3 

0. 

3. 

3. 

3. 

3. 

3. 




totals 

8.507E-38 

7. 850E-08 

9.228E+01 

3. 

3. 

3. 





DETECTOR NJM3ER 
COORDINATES 

1 

X = 3. 

V 

= 0 . 

Z = l. 

353335+01 

L34ER 

ENERGY = 2.003E-01 ME V 


UPPER ENEGY 
BDSIMEV) 

UFLUX 

ST3 DEV 

O DEV 

SF LUX 

STO 05V 

PC OEV 

TFLJX 

STD DEV 

PC OEV 

4. OOOE-31 

3. 

0 . 

3. 

6.967E-05 

L. 2335-35 

1.342E+01 

6.967E-05 

1.283' -05 

1.84iE+3l 

6. 030E-31 

3. 

0 . 

3. 

7. 7425-35 

2.1935-35 

2.327E+01 

7.742E-05 

2 . 18ot-0 5 

2.3276+01 

1 .OOOE+OD 

3. 

3. 

3. 

3. 796E -05 

6. 3 695 -36 

1.310E+01 

3. 7 96E-05 

6.8591-06 

1.81 > E + 3 1 

1 . 500E + 3 3 

0 . 

0 . 

3. 

3.476E-05 

5.3535-36 

1.540E+01 

3.476E-05 

5.3536-06 

1.54 ■ 6 + 3 L 

2 • 000E+03 

3. 

0 . 

3. 

3.323E-35 

2.9725-36 

9. 315E+00 

3.328E-05 

2.972E-06 

9.31^6+00 

2-500E+30 

3. 

3. 

0 . 

3. 43B5-35 

5.6295-35 

1.537E+31 

3.43BE-05 

5.6?9r- 06 

1. 637E+01 

2 • 999E+03 

0. 

0 . 

0 . 

2.935E-35 

3. 3365-35 

1. L25E+01 

2.935E-05 

3.306E-06 

1 .126E+01 

3.0D0E+33 

2-649E-34 

2 .98 BE-35 

1. 128E + 01 

3.445E-33 

1.0985-03 

3. 136E+01 

2.650E-04 

2.9P8f-05 

1.1286+01 

TOXALS 

2.649E-34 

2.988E-35 

l. 128E+01 

3.1395-34 

2.6471-05 

8. 435E+00 

5.783E-34 

3.9925-05 

6. 393E + 00 
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GAMMA DOSE 

RATE IR/HR 

PER PH3T3N/CM2 

-SEC) 





UPPER ENEGY 
8DS{ ME V ) 


STD DEV 

' PC DEV 


ST3 DEV 

PC 3EV 


STD DbV 

3 C DEV 

4.3D0E-31 

D • 

o. 

3. 

4.222E-U 

7. 7775-12 

U34ZE+31 

4.222E-11 

7.777u-12 

1. 34'E+31 

&.ODOE-31 

0. 

3. 

3. 

8.149E-U 

2.3335-11 

2- 3276+01 

8. 149E-11 

2.33 lb-11 

2.8275*31 

1.000E+33 

0. 

0. 

3. 

6.1225-11 

1.1385-11 

1.313E+31 

6.122E-11 

l.llMc-ll 

1.81 . f *31 

l • 500E-+ 0 3 

3. 

o. 

0. 

7.931 E- 11 

l. 2175-11 

1 . 5 4 3E+ 31 

7.901E-11 

1.21 7b- 1 1 

1.54 t * 3 L 

2. 0006+03 

3. 

0. 

3. 

8.7775-11 

B. 5145-12 

9. 31 5E + 3 3 

8.777E-11 

8.614t-12 

9. 81 1l*3D 

Z.500E+33 

0. 

0. 

o. 

1.2DSE-1D 

1.9755-11 

1 • 537E + 31 

1.206E-13 

1 .9 7*^E-1 1 

1.637 b*3l 

Z.999E+D3 

0. 

0. 

3. 

1.198E-13 

1.3535-11 

J.125E+01 

1.19BE-13 

1.35uE-ll 

1.12 t*3l 

3.300E+33 

1.162E-39 

1.311E-13 

1.L28E+31 

1.5115-13 

4.8155-14 

3. 1 3 5E + 31 

1.162E-39 

1.3 lib- 13 

1.1235*31 

TOTALS 

1.162E-39 

1.311E-10 

1.128E+31 

5.923E-13 

4.9965-11 

8 • 435 E + 3 3 

1.754E-39 

1.2L3E-10 

5.3986*33 


DETECTOR NJM3ER 
COORDINATES 

2 

X = 0. 

Y 

* 0 . 

Z = 2. 

13333 E* 31 

-OVER 

E M E R GY = 2.333E-01 MEV 


UPPER ENEGV 
BDS(MEV) 

UF LUX 

STD DEV 

»C DEV 

sflux 

SO D5V 

PC DEV 

TFLUX 

STD DLV 

PC DEV 

4.000E-DI 

0. 

0 . 

3. 

1. 4786-35 

3. 7585-35 

2.542E*3i 

1.47BE-35 

3.75Bb-36 

2. 54 'f +31 

6.000E-31 

3. 

3. 

0 . 

5.7335-35 

1.3595-35 

1.375E*31 

5.703E-36 

1.369E-06 

l. 875E+31 

l.OOOE+33 

0 . 

0. 

3 • 

1.572E-05 

3.1545-36 

5.1 37E*0l 

1. 572E-35 

8.154F-06 

5.1875+31 

1.5006*03 

0- 

0 . 

0 . 

4.839E-36 

8.4165-37 

1.739E*3l 

4. 839E-35 

8.41 t.F-37 

1. 7396+31 

2.0306*33 

3. 

0 . 

0 . 

6- 2835-35 

1.9445-36 

3.395E*31 

6.283E-35 

1 .9446-06 

3. 3956 + 3L 

2.500E*03 

0 . 

0 . 

o. 

5.4975-35 

7.3385-37 

1.421E*31 

5. 497E-05 

7 .8 Do E~0 7 

1.4216*31 

2 .999E+03 

0 . 

0 . 

3. 

5.86BE-05 

5.3545-07 

1.3836+3 l 

5.B6BE-35 

6.354L-37 

1.383E+31 

3. 300E*33 

2.4326-35 

2.7436-06 

1.128E*01 

6.1936-39 

1.9335-39 

3.373E+01 

2.433E-35 

2.74 1C - 06 

1.1216+31 

TOTALS 

2 . 432 E-3 5 

2.743E-06 

1.128E*31 

5.3735-35 

1.3135-35 

1.721E+31 

8.302E-35 

1 .347 C - 05 

1.261E+31 




GAMMA DOSE 

RATE tR/HR 

PER P H0T3M/ C M2 

-SEC) 





UPPER ENEGY 
BDSIMEV) 


STD DEV 

PC DEV 


SO DEV 

PC DEV 


STD DEV 

PC DEV 

4.0006-01 

3. 

0 . 

0 . 

B.959E-12 

2.2735-12 

2.542E+31 

8.959E-12 

2.278 — 12 

2. 54^ E + 31 

6.000E-01 

0 . 

0 . 

3. 

6.3005-12 

1.1255-12 

1. 375E+01 

6.303E-12 

1,12.6-12 

1.875E+31 

1.000E+03 

0 . 

0 . 

0 . 

2.535E-11 

1.3155-11 

5. 137E+01 

2. 535E-1 1 

1.310E-U 

5.1875+31 

1.5006*33 

3. 

0 . 

0 . 

1.1335-11 

1.9135-12 

l. 739E+31 

I.IOOE-H 

1.913b- L2 

L. 7336+31 

2 • OOOE+-03 

3. 

0- 

3. 

1.821E-11 

5.5365-12 

3.395E+DI 

1.8 21E-1 1 

5.636E-12 

3. D95E+31 

2.500E+03 

0 . 

0 . 

3. 

1.9295-11 

2.7405-12 

1* 421E + 31 

1.929E-11 

2.74 F-12 

1.4216+31 

2.999E+33 

0 . 

0 . 

o. 

2.395E-1L 

2.5935-12 

1. 333E+D1 

2.395E-11 

2.593E-12 

1. 3836+31 

3.300E+33 

1.367E-10 

1.203E-11 

1.123E+31 

2.7165-14 

3.3475-15 

3. 373E+31 

1.367E-13 

1.203.-11 

1. I23t+3l 

TOTALS 

1.367E-13 

1.203E-11 

1.1285+31 

1.128E-13 

1.3415-11 

1. 721E+31 

2.195E-13 

2.76 7fc- 1 1 

1.26U + 31 
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DETECTOR NJM3ER 
COORDINATES 

3 

X = 0. 

T 

= 0. 

Z = A. 

233336*31 

LOWER 

ENERGY = 2.333F-01 

MEV 


UPPER ENEGV 
BDSIMEVJ 

jflux 

STD DEV 

>C DEV 

SFLUX 

ST3 3; V 

p: 3EV 

tfljx 

STD DEV 

DEV 

A.030E-31 

0. 

0. 

3. 

1.97BE-3S 

1.5775-35 

7. 972E+31 

1.97BE-35 

1 .577 f -36 

7.9726*31 

6.000E-01 

0. 

0. 

3. 

2.391="07 

9. 9355-38 

A. 155E*3 l 

2.391E-37 

9 .9 3 > 

i -38 

A. 1556*31 

1 • OOOE *00 

0. 

0. 

3. 

A.78BE-07 

1.5735-07 

3. 23AE01 

A.78BE-37 

1.573 

-3 7 

i.28Vt*3l 

1 . 500E + )3 

0. 

o. 

3 • 

2.1A7E-37 

3.9525-33 

1. 3 A1E * 3 l 

2.1 A7E-37 

3.95 

£-3 8 

1 . 3 A 1 * f 3 1 

2. 030E *0 3 

o. 

3. 

3. 

1. B3BE-37 

3. 3AAE-39 

2 . L27EO 1 

1.803E-37 

3 . 8 A ■: 

t - 38 

2. 127031 

2 . 500E* 33 

3. 

3. 

3. 

5.772E-07 

2.3225-37 

A. 323E*31 

5.772E-37 

2.322 

L -0 f 

A. 32 *6*31 

2.9996*33 

0. 

0. 

3. 

A.A58E-37 

1.3525-07 

3. 355EO! 

A. A59E-37 

1.36 t 

L-3 7 

3.356c *31 

3.000E+00 

B.2D3E-37 

9.2A9E-3B 

1.128EOI 

5.8735-13 

3. 4135-1 3 

A. 96lE*3l 

8.207E-37 

9.2A^, 

l- 38 

I. 127i *31 

TOTALS 

8.230E-37 

9.2A9E-0B 

1.128EOI 

4. 1 1 5E -3 5 

1.5555-35 

3. 732EO! 

A.935E-35 

1 .559 

-06 

3. 15 ^E*31 




GAMMA DOSE 

RATE U/HR 

PER PHOTON/ CM2 

-SEC) 





UPPER ENEGY 

bdsimevi 


ST3 3EV 

’C 3EV 


ST3 DEV 

=>; DEV 


STD DEV 

=>C DEV 

a.OOOE-31 

0. 

3. 

3. 

1. 1996-12 

9.555E-13 

7.972E*31 

1.199E-12 

9.556T-13 

7. 972E*31 

6. 0D0E-31 

0. 

0. 

D. 

2.5155-13 

1 . 3 A5E-13 

A.155E*31 

2.516E-13 

1 .0A6L-13 

A. 156E*31 

1 . 030E* 03 

0. 

0. 

3. 

7.7225-13 

2.535E-13 

3.234E*01 

7.7 22E-13 

2.536l-13 

3. 29 At* 3 L 

1.500E+33 

3. 

3. 

0. 

A. B795-L3 

3.9325-1A 

1.3A1E*31 

A.879E-13 

8.982014 

1. SAlt*31 

2.000E*0 3 

0. 

0. 

3. 

5.239E-13 

1.11 AE-13 

2. 127E*31 

5.239E-13 

1 .114013 

2. 1 27 E*3 1 

2 • 500E* 33 

o. 

o. 

3. 

2.325E-12 

8. 149= -13 

A. 323E*D1 

2. 3 25E-12 

8 . 1 A9 r - 1 3 

A. 3236*31 

2 • 999E+3D 

3. 

0. 

3. 

1.82DE-12 

5.550E-L3 

3. 35 5E * 3 1 

1.823E-12 

5.56 013 

3. 355E * 3 L 

3 . 0 00E*03 

3.597E-12 

A.D57E-13 

1.128c *31 

3.31AE-15 

l. A95E-15 

A.951E*3L 

3.503E-12 

A.057E-13 

1.127EOI 

totals 

3 . 59 7E-1 2 

4.D57E-13 

1. 1285*01 

7.0B2E-12 

2. 579E-12 

3. 732 E*3 1 

1.3 68E-11 

3.37At-12 

3.159E*31 


END OF PROB.EM 20 

FINAL RAN031 NUMBER 
056055355315 

*01* EXIT IN EXECT 
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Figure 1. - Geometry illustration. 
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Figure 2. - COHORT-II major subdivisions and functions. 


62 

















Figure 5. - Subroutine EXECT. Acts as executive routine for entire code. 


IHI = NGY 



Figure 6. - Subroutine XSEC. Given the particle energy (EA) P performs linear interpolation 
to obtain a cross-section value at that energy. 


















Figure 7. - Subroutine MAIN. Computes storage requirements of problem and reads in part of input data. 

































Figure 9. * Subroutine SORCE. Main routine for primary particle source ge- eratcr. 
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Begin loop on 
spatial coordinate 




Figure 10. - Subroutine COORD. Selects spatial coordinates of primary source particles according to input spatial distributions. 
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AZWT = 1. 0 


Compute weight 
from bias *- 
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Compute overall 
angular weight 
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ANGWT = PWT*AZWT 


( /stop ) 


Figure 11. - Subroutine DIRCOS. Generates direction cosine for source particles according to option. 




















re 12. - Subroutine ENERGY. Selects source particle energy from energy spectrum. 
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Figure 13. - Subroutine HTA. Executive routine for particle tracking and scoring part of the code. 
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Figure 14. - Subroutine HISTRC. Performs particle tracking through reactor shield geometry. 
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Figure 15. - Subroutine SCATR. Performs particle collision mechanics and calculates scattered flux at point detectors. 
















Figure 16. * Subroutine GEOMT. Performs geometry calculations. 


Figure 17. - Subroutine OUTPUT. Writes out tables for energy deposition in regions. 













Figure 18. - Subroutine TLFLX. Normalizes track length fluxes and writes leakage tape. 



Figure 19. - Subroutine FLPRT. Writes out region fluxes and particle currents leaking from shield geometry. 


















Figure 20. - Subroutine F PRINT. Prints out point detector fluxes. 













Figure 2L - Subroutine SECGAM. Computes and allocates storage required for secondary gamma problem. 



Figure 22. - Subroutine SECRED. Reads in secondary gamma input data. 
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Figure 23. - Subroutine SECSO. Computes secondary gamma and writes source tape. 
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